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AN APPLICATION OF CROCCO’S STREAM FUNCTION 
TO THE STUDY OF ROTATIONAL SUPERSONIC 
FLOW PAST AIRFOILSt 


By ABRAHAM KOGANtY 
(Technion—Israel Institute of Technology, Haifa, Israel) 


[Received 25 October 1956] 


SUMMARY 

A method of successive approximations to the ideal flow around airfoils at high 
supersonic Mach numbers is developed, based on the concept of Crocco’s stream 
function. 

Taking the undisturbed flow as a starting-point, a second approximation to 
the flow in the whole region surrounding the airfoil is derived, and the results for 
velocity and pressure distribution at the airfoil surface are compared with the 
corresponding expressions obtained by potential flow theory. 

The flow behind a plane shock wave is next chosen as zero-order approximation 
in order to obtain better results for the flow in the vicinity of the leading edge. 
The first approximation gives expressions for the shock-wave curvature and for the 
derivative of the pressure coefficient at the leading edge, which check similar results 
obtained by other investigators. The second approximation yields expressions for 
the derivative of the shock-wave curvature and for the second derivative of the 
pressure coefficient at the leading edge. 


Nomenclature 
A, B,C constants defined by equations (78). 
G,H,L constants defined by equations (56). 


a;;, 54; coefficients in expansions of velocity perturbations. 

c speed of sound, in units of limiting velocity. 

C, specific heat at constant volume. 

Cy pressure coefficient. 

g(a) airfoil profile. 

G(x) airfoil profile form function. 

K,,,, Kg, ®irfoil and shock curvature at leading edge, respectively. 
M Mach number. , 

P function of w, defined by equation (11). 

P;; coefficients in expansion of P in powers of , and y,. 


+ Abstracted from a doctoral dissertation submitted by the author in the Department 
of Aeronautical Engineering, Princeton University. The work was sponsored by the U.S. 
Air Force, under contract AF61(514)—870. 

The author gratefully acknowledges the helpful suggestions of Professor Sydney Goldstein, 
under whose supervision this research was carried out. 

t Senior Lecturer, Department of Aeronautics. 
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pressure. 
entropy. 
temperature. 
cartesian coordinates. 
components of w in x- and y-directions. 
components of velocity perturbations. 
velocity, in units of limiting velocity. 
constants defined by equations (81). 
\(Mj—1). 
« (M3—1). 

, A,, A, determinants defined by equations (88). 
flow inclination to undisturbed stream. 
thickness parameter. 

Mj sin’. 
coefficients in expansion of 7 in powers of ¢ or ‘V. 
sin-1(1/M). 
shock-wave inclination to undisturbed stream. 
density. 
length parameter along shock profile. 
slope of shock-wave profile. 
Ty: Ty, T». Coefficients in expansion of r in powers of x. 
a Crocco’s stream function. 
ys perturbation stream function. 
ys, wy,... terms in expansion of yw. 
Vo, do shock and flow inclination at leading edge. 
My, Ug -—« zero-order approximations to M and uw. 
Po, So. Ty, Pp stagnation values of p, S, 7’, and p. 
Py, Sy. Ty, py, My, Uy, 01, Wy, ¢, values of respective parameters in undis- 
turbed flow. 


Po, So, To, po, Mg, Ug, Vg, Cg Values of respective parameters behind plane 
shock wave. 


1. Introduction 

THE potential flow approximation to the supersonic flow around airfoils 
gives a fair estimate of the actual flow when the airfoils considered are thin 
and the free stream Mach number small enough. However, when the thick- 
ness ratio of the airfoils and the free stream Mach number are increased, 
the vorticity generated by the curved shock wave attached to the leading 
edge becomes an important factor, and approximations obtained by poten- 
tial flow theory become rather poor. In fact, when the thickness ratio and 
the free stream Mach number are increased to the point at which the free 
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stream Mach angle equals the tip wedge angle, potential flow perturbation 
solutions lose their physical meaning altogether (1). 

Since the method of characteristics, which is capable of giving exact 
solutions of the ideal flow, is rather lengthy and provides results only in 
numerical form, attempts have been made in recent years to devise approxi- 
mate methods which should simplify the treatment of high supersonic flows 
without neglecting the main rotationality effect (2—5; see also 16 and 17). 

The study of flow behind a curved shock wave in the vicinity of the 
leading edge of airfoils is particularly interesting, since essential features 
of the high supersonic flow may be explained by such studies. Also the 
approximate methods for the calculation of the whole flow field take the 
flow in the region of the leading edge as a starting-point. 

The aim of the present investigation is to develop a method of successive 
approximations for the study of flow at high supersonic Mach numbers, 
based on the use of Crocco’s stream function. 

A first and second approximation to the flow in the whole region sur- 
rounding the airfoil are calculated, starting from the undisturbed flow as 
the zero-order approximation. The well-known results obtained by Buse- 
mann (6, 7) and Van Dyke (1) by potential flow theory are obtained, but 
it is shown that their validity extends to much higher values of the free 
stream Mach number and of the leading edge angle. 

Then taking the flow behind a plane shock wave as zero-order approxima- 
tion, the flow in the vicinity of the leading edge is studied. The solution 
of the first approximation gives expressions for the shock-wave curvature 
and for the slope of the pressure coefficient at the leading edge which check 
the results obtained by the method of characteristics (8-12). The second 
approximation yields formulae for the derivative of the shock-wave curva- 


ture and for the second derivative of the pressure coefficient at the leading 
edge. 


2. Crocco’s equation for two-dimensional flow 
Consider a sharp-nosed airfoil placed in a uniform supersonic flow of an 


ideal and perfect gas. Crocco (13) showed that it is possible to define a 
stream function ‘Y by the equations 


¥, ots u(1—w?)iy-D 
, Nap. 
where suffixes denote derivatives. 


The continuity equation is satisfied identically by ‘Y, and the momentum 
and energy equations lead to Crocco’s equation 


(c?@—u?)¥,,— 2u0'¥,, + (c2—v2)'F,,, = (1—w®)'-D(w2_ 2) f(W) (2) 


(1) 
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p pressure. 

S entropy. 

T temperature. 

ry cartesian coordinates. 

u,v components of w in x- and y-directions. 

u',v’ components of velocity perturbations. 

w velocity, in units of limiting velocity. 

Xy, %, %3 constants defined by equations (81). 

x a (M}—1). 

B / (M}—1). 

A,, Ay, As, Ay, A; determinants defined by equations (88). 
8 flow inclination to undisturbed stream. 

€ thickness parameter. 

n Mj sin?v. 

No» Ny». coefficients in expansion of » in powers of ¢ or ¥’. 
My sin-1(1/M). 

v shock-wave inclination to undisturbed stream. 
p density. 

o length parameter along shock profile. 

T slope of shock-wave profile. 

To: Ty, Tg Coefficients in expansion of r in powers of x. 
¥ Crocco’s stream function. 

M7 perturbation stream function. 

y,, we,... terms in expansion of w. 

Vo, do shock and flow inclination at leading edge. 

My, Uo zero-order approximations to M and u. 


Po. So, To, Po stagnation values of p, S, 7’, and p. 
Py, Sy, Ty, py, My, Uy, %, Wy, ¢, values of respective parameters in undis- 
turbed flow. 


Pe, So, To, po, Mz, Ug, Vg, Cg Values of respective parameters behind plane 
shock wave. 


1. Introduction 


Tue potential flow approximation to the supersonic flow around airfoils 
gives a fair estimate of the actual flow when the airfoils considered are thin 
and the free stream Mach number small enough. However, when the thick- 
ness ratio of the airfoils and the free stream Mach number are increased, 
the vorticity generated by the curved shock wave attached to the leading 
edge becomes an important factor, and approximation: obtained by poten- 
tial flow theory become rather poor. In fact, when the thickness ratio and 
the free stream Mach number are increased to the point at which the free 
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stream Mach angle equals the tip wedge angle, potential flow perturbation 
solutions lose their physical meaning altogether (1). 

Since the method of characteristics, which is capable of giving exact 
solutions of the ideal flow, is rather lengthy and provides results only in 
numerical form, attempts have been made in recent years to devise approxi- 
mate methods which should simplify the treatment of high supersonic flows 
without neglecting the main rotationality effect (2-5; see also 16 and 17). 

The study of flow behind a curved shock wave in the vicinity of the 
leading edge of airfoils is particularly interesting, since essential features 
of the high supersonic flow may be explained by such studies. Also the 
approximate methods for the calculation of the whole flow field take the 
flow in the region of the leading edge as a starting-point. 

The aim of the present investigation is to develop a method of successive 
approximations for the study of flow at high supersonic Mach numbers, 
based on the use of Crocco’s stream function. 

A first and second approximation to the flow in the whole region sur- 
rounding the airfoil are calculated, starting from the undisturbed flow as 
the zero-order approximation. The well-known results obtained by Buse- 
mann (6, 7) and Van Dyke (1) by potential flow theory are obtained, but 
it is shown that their validity extends to much higher values of the free 
stream Mach number and of the leading edge angle. 

Then taking the flow behind a plane shock wave as zero-order approxima- 
tion, the flow in the vicinity of the leading edge is studied. The solution 
of the first approximation gives expressions for the shock-wave curvature 
and for the slope of the pressure coefficient at the leading edge which check 
the results obtained by the method of characteristics (8-12). The second 
approximation yields formulae for the derivative of the shock-wave curva- 


ture and for the second derivative of the pressure coefficient at the leading 
edge. 


2. Crocco’s equation for two-dimensional flow 
Consider a sharp-nosed airfoil placed in a uniform supersonic flow of an 
ideal and perfect gas. Crocco (13) showed that it is possible to define a 
stream function ‘VY by the equations 
x ; (1) 
WY, = —v(1—w*)o—) 
where suffixes denote derivatives. 
The continuity equation is satisfied identically by ‘Y, and the momentum 
and energy equations lead to Crocco’s equation 


(c?—u®) ¥,,— 2uv'¥, +(8—v7)¥, = (L—w?)'tO-D(wt—c?) f(¥) (2) 
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, 1 dS ‘ 
where f(¥) = yc, de" (3) 

The only restrictions made in deriving (2) are the assumptions of a two- 
dimensional steady flow and of a frictionless perfect gas. Equation (2) is 
valid for compressible and rotational flow. It is suitable therefore for the 
treatment of two-dimensional flow behind shock waves. 

It is remarkable that (1) define ‘Y only by velocity components. No 
thermodynamical state property of flow appears explicitly in them. This 
feature of Crocco’s stream function gives it the advantage over the ordinary 
stream function for compressible flow. It is also interesting to observe that 
since by the energy equation 


” 
SiN an ea (4) 


equations (1) may be written in the form 


W W(y—-)) —S, 

= (= GP 
T\'y-D S—S, 

reel SE 
T) ‘ R Po 


so in the case of isentropic flow Crocco’s stream function becomes identical 
with the ordinary stream function for compressible flow. They are really 
different only in the case of a rotational flow. 

Equations (1) may be written in the form 





d¥ = (1—v*)"-Y(udy — vdz), (5) 


which shows that ¥ is a constant along streamlines and that it increases 
continually when we move from the airfoil surface towards more remote 
streamlines; ‘Y can therefore be used as a coordinate along the shock wave 
profile. In particular, the shock-wave parameter 7 = Mj sin*v is a function 
of ‘VY: . . 
7 = 7(¥). (6) 
It is convenient to express f (‘’) in (3) explicitly through ». The following 
equations relate the entropy increase through a shock wave to the shock- 
wave parameter 7: 
s—2, 


v 


p 2y ( ean A, ¥— ot 25) />. 
— = —_|n———]}; — * 
Py yt! 2y p yt+l , 





= In(p/p,)+y In(p,/p), 
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Using these relations, we obtain for f(‘¥) 


at See 
I) = 5 P= ae (7) 





with Qn) = - (8) 
aln—{(y—1)/2y} n+ {2/(y—1)}] 


3. The perturbation stream function 
Define a perturbation stream function ¥(z, y) by 


YF = u(1—up)!-fy+-Y(a, y)}, (9) 
where uy is a constant representing the zero-order approximation to wu. 
Expressed in terms of #, (2) becomes 


€ 


= i (i— ui) —{(c?@—u* yp, — 2urp py + (c?—v* ph, y} = PwIQn) F, ; 
(10) 
where P(w) = — Bia ug 1(1— ud) -D(1 — w? 7 + Dy -D(w?—c?), (11) 
vly—1) 


We shall choose the coordinate system with the z-axis in the direction 
of the uniform flow which represents the zero-order approximation. The 
perturbations of the velocity components u’ = u—w, and v’ = v can be 
expanded by (1) into power series of the first derivatives of ¥ 


Wy = > ay Vid (12) 
v' |g = Par by oh. (13) 


The coefficients a,; and b;; may be expressed through the Mach number M, 
of the zero-order flow approximation. The coefficients up to the second 
order are 

Log = Igy = Ay = bo = bq = Faq = Fog = 9) 

Ayo = —1/(Mp—1) 

og = {MG/(ME—1)8}{(2—y) ME—3} 

Gog = —3MG/(MG—1) 

b = —1 

by = MG/(MS—1) ] 


Equations (12) and (13) will be used subsequently in order to expand 
the coefficients of (10). 


(14) 
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4. Boundary conditions 
If the airfoil profile is given by 


y = g(x) (15) 
the condition of tangency of flow at the solid boundary may be expressed 
by ¥(z, g(x) = 0. 


Or, taking the derivative of this equation along the airfoil profile and 
introducing the perturbation stream function y, the tangency condition 


— Ho(x, 9(x)) +921 +Hy(x, g(x))} = 0. (16) 
The Rankine-Hugoniot shock-wave conditions may be expressed as 


follows: W,—Wy = 9, (17) 
9 


yi tin) (18) 


Wi n(W,— Wy) = 


where w, and w,, are the components of the non-dimensional velocity in a 
direction tangential and normal to the shock-wave profile respectively, 
evaluated right behind the shock wave, and the quantities marked by the 
index , are evaluated in front of it. 

In terms of the velocity components wu, v in the x, y coordinate system, 
and the slope 7 of the shock wave, w, and w,, can be expressed as 





U-+-TV 
wo, = Tee Ss (19) 
‘ f(i+78) 
TuU—v 
— é (20) 
\ (1+-7?) 
When (19) and (20) are inserted in (17) and (18) these become 
utry = u,+7r, (21) 
9 ~~ 
and (ru, —v,)(ru—v) = ——c?(1 4-28) 4 Yo! re —v,)*. (22) 
y+! y+1 


Our problem is to find a solution (x, y) of (10) subject to the conditions 
(16), (21), and (22). The differential equation (10) is of the second order 
and a solution would be determined by two boundary conditions. Yet the 
conditions (16), (21), and (22) do not over-determine the system since they 
contain a new unknown, the shock-wave parameter » = 7(‘’) or tr = 7(z). 
The problem is further complicated by the fact that the right-hand side 
of (10) depends itself on the unknown parameter 7. 

The natural approach seems to be one of successive approximations. 
We shall assume a zero-order approximation for the velocity field and for 
the shock profile. Using this approximation, we shall express explicitly 
the differential equation and the boundary conditions correct to the first 
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order. This will determine a first-order approximation to the flow field 
and the shock-wave profile. Then the procedure will be iterated in the 
hope that we thus obtain a sequence of solutions which converge to the 
correct result. 
Up to this point we have assumed that the zero-order approximation to 

the flow field is given by a uniform flow 

W = Uy 
and consequently by a plane shock wave 

T=T) 
without specifying, however, the values of wu, and 7». We can expect the 
results of higher approximations to be better the closer the zero-order 
approximation itself represents the flow. In the case of very thin airfoils 
we may hope to obtain good results by taking the undisturbed flow in front 
of the airfoil as the zero-order approximation. But if the airfoil is thick, 
the free stream would provide a bad starting-point. In this case the results 
will be much better if we choose as zero-order approximation the flow 
behind a straight shock-wave generated by a plane wedge of wedge angle 
equa. to the leading edge angle of the airfoil. We shall consider separately 
these two cases. 


5. First approach: free stream taken as zero-order approximation 
According to our convention above, we choose in this case the z-axis in 
the direction of the undisturbed flow. 
Let us introduce a thickness parameter « by writing the airfoil profile 
equation in the form 


y = g(x) = € G(x), (23) 
where O(2) nex = 1. 
The shock-wave parameters may be expanded in powers of e: 
» = M?sin*vy = yy +7, +e +... (24) 
or, alternatively, 7 = tanv = 7)+7,€+7,€7+.... (25) 


The perturbation stream function may be similarly expanded 
f= ppetyee+.... (26) 


In the present case the zero-order approximation to the second shock 
condition (22) reduces to 


1+73 
M} = —_* 
2 > 
or, defining vg by tanvy = Tp, 
Pace. 
1 sin®y, 


This yields No = Mj sin*y, = 1. (27) 
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But since dS/(c,dn) contains the factor (y7—1)? and dy/d¥ = O(e), it 
follows that the right-hand side of the differential equation (10) is at least 
of third order in e. The first- and second-order expressions for the differen- 
tial equation reduce therefore to 


Pryy— Orr = 0 (28) 
and 


sed 


a0, 1 
Pry Pr22— 2(a3+ I),, Pizy— 


iT 
—(/+1)— ma 


Poyy— pore = —(y+1) 


respectively, where a? = Mj—1. 
The tangency condition (16) will be expanded around y = 0. Inserting 
(23) in (16) and equating equal powers of «, we obtain for the first approxi- 


mation wh, -(2, 0)+G"(a) = 0 (30) 
and for the second approximation 
Po2(x, 0)+ G(x)by 22, 0) - G' (x)byy(x, 0) = 0. (31) 


The shock-wave conditions should be expressed at the shock-wave profile. - 
But since this is unknown as yet, we may expand quantities at the actual 
shock line in Taylor series around the zero-order position of the shock line. 
We thus obtain from (12) and (13) the velocity components on the shock- 
wave surface as 

Zz 


ule, r dr] 


0 


= af} +4yo%yy€4 la | 74(x) dx ayy ro Pay-+ Aan by + doa el to] 


0 
(32) 
and 


of, [ T as] 


0 


\s a: 
-— 2 o1¥i2€+ box f 7 (x) dx Vary tbo Port Ons Piz Pay + [> (33) 
0 
where y,,, iy, Yizy and y,,,, are evaluated at y = t)x. Since in the present 
case % = %, vy = 0 
the shock-wave conditions (21) and (22) simplify to 
u v 


Suet (34) 
Uy Uy 
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u v 2 @& 2 —_ 
and 7 — Y i} 35 
5S ~ ytluj wef gees asi 
Inserting here the expansions (32), (33), and (25), we find that (34) is 
identically satisfied to the zero-order of approximation, and that it gives 


to the first orde 
. Ay yy +79 Oo, 4, = 0 (36) 





and to the second order 
z 
B49 Poy +7 O01 Port (Ayo PryytTo bo Przy) | 74(x) dx + 
0 


+bu 71 Pir +9 Pi y+Ooe Pir +794; Pir Pry =90. (37) 


Equation (35) gives to the zero order 


‘aa 
1 sin2y, 
to the first order 
To 0 Piy— on ti: = aald- 1}n (38) 
y+l\uy 


and to the second order 
2 c 
y+ (1 = G)tet+ 21072) + 18 tio bey Foon Port 
+(7§ Ay Piyy—To bo Piry) | 7,(x) dx + 
0 


+( 279490 $1 —O01 Piz)t1 +75 G20 Ply +76 202 Fiz — 
—T boy tiety = 9. (39) 
Here all derivatives of % are evaluated at y = 792. 


The first-order problem is defined by the differential equation (28) with 
the boundary conditions (30), (36), and (38). The general solution of (28) is 


f(x, y) = m,(x+ay)-+n,(x—ay), 
where m and n are determined by (30) and (36). We obtain 


vy(z,y) = —G(x—ay). (40) 
Using this result in (38), we obtain for 7, 


y+1 (a? - @'(0) 


ep ll 41 
ee) Oe (41) 
Hence to the first order the shock wave is still plane, but it no longer 
coincides with the Mach wave. 
The second-order problem is defined by the differential equation (29) with 
the boundary conditions (31), (37), and (39). These equations are deter- 


mined by inserting for %, and 7, the first-order solution (40) and (41). 
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Equation (29) becomes 


= tte “(a —ay)G" (x —ay). (42) 


Poyy—@Worr -_ (y+1 ) 
The following particular integral J, of this equation was found by 
inspection: 41 (a2-+-1)2 
.= $e (a > ) {G’(a—ay)}*x. 
4 a 
The general solution of (42) may be expressed, therefore, in the form 


yt} (f+ IP gq 
3 


Pale, y) = mylx-+ay)+n4(2—ay) — "= (x—ay)}2a. 


The application of the shock-wave condition (37) yields 
m,(x) = 0, 
and we can put, without loss of generality, 
mM, = 0; 
n, is then determined by the tangency condition (31) and we obtain 
_ yt] (o® +1) 


f(x,y) = 4 72 


{G'(x—ay)}?y—a G(x—ay) G'(x—ay). 
(43) 
Finally, inserting this value of , into the second shock-wave condition 
(39), the follewing expression for 7, is found: 


mt 2 3y—1 
7,(x) = (y os (a?- +1) ( 42 9 By 8 +5)(@"(0))*+ 
32 *: y+1 


+5. 480(0)"(0)-+2(a2-+1°6'(0) 6 “(oy (44) 
7 


To the second ieee the shock wave is represented by a parabola. 
Now using the expressions (40) and (43) for %, and #,, we can calculate 
u/u, and v/u, at a point x, y of the flow, correct to the second order. They 
depend in general on G(x—ay), G’(x—ay), and G"(x—ay). But at the 
airfoil surface the terms of u/u, containing @” cancel and we obtain 

u(x, ¢ Me)! G'(e)e— 71 a+ 2y— Wy + Wia® +1 sar gyrags 


4 
Uy a a 





(45) 
This result is identical with the one obtained by Van Dyke (1) by potential 
flow theory. 
Inserting this and the corresponding expression for v/u, into the expanded 
expression of c,,, we obtain to the second order 


Cy = 2 (x)e 4 LEMME TI 40 Gr (ayi2e, (46) 
a aM 
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which coincides with the well-known result of Busemann (6, 7). Taking 
the free stream as zero-order approximation, the pressure coefficient at the 
airfoil surface depends, to the second order, only on the local slope. 

As far as conditions at the airfoil surface are concerned the present 
method does not improve on the results of potential theory. This is due 
to the fact that no terms representing vorticity appear in the differential 
equation up to the second-order approximation, and the boundary con- 
ditions at the airfoil are also equally approximated by both methods. 

But our method also gives good results away from the airfoil, even near 
the shock wave, where the potential theory fails because of its inherently 
incorrect application of the shock-wave conditions. A check on our results 
(27), (41), and (44) for the shock profile may be obtained by expanding the 
exact shock profile equation of the flow over a plane wedge of wedge angle 





5 = tan-'e. 
The exact shock profile equation in this case is 
§ al EE M; ti |; 
€ 2 {M2 7?/(1+7?)—1} , 


where 7 is the tangent of the shock angle. Inserting 
T = Tot7,€+7,€7+... 
and equating equal powers of e, we obtain 








oii =, (47) 
ee es (48) 

_ (y+? (a? +1)*f Sy—1 ., gl 
ae eT eee wa 


Equation (47) corresponds with (27). When we insert (47)—(49) in (41) and 

(44) we find that we must have 
G’(0) = 1 
and G’"(0) = 0. 

In this connexion it is interesting to observe that in potential flow theory 
the region of disturbed flow is bounded, to any order of approximation, by 
the airfoil surface and by the Mach wave from the tip in the undisturbed 
flow. And when M, and the*tip angle 8, are so high that 


1 
V(MZ-1)’ 
the region of definition of disturbed flow shrinks to zero and potential flow 
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approximations become altogether meaningless. Such a breakdown occurs, 
to be sure, in the present approach too, but for a given tip angle it is delayed 
to much higher values of M,. Thus in the present theory the ratio of the 
tangents of the shock and airfoil tip angles is given to the first approximation 
by 4 

. l 4ytl_ Mt (50) 


tand, ¢«@(0)/(M?—1). 4 (M?—1)’ 
and in potential theory by 





T ] 


tand,  «@’(0)(M?—1) (51) 

















TABLE | 
T 
M which =I 
1 for whic and, 
8 Potential theory | Present theory 
_ 11°48 28-8 
10° 5°75 14°0 
ss 3°86 g:08 
20° 2°92 6°58 
25° 2°36 4°75 











In Table | are indicated the values of M, for which these ratios become equal 
to 1 and breakdown of each theory occurs. It shows that Busemann’s result 
for c,(x) at the airfoil surface is of value beyond the region for which it was 
originally calculated. 

There are no added difficulties in the calculation of a third-order approxi- 
mation than in the preceding calculation, except for some more lengthv 
algebra. The third-order approximation to the differential equation will 
contain a term arising from the right-hand side of (10). This term will 
depend only on the first-order quantities %, and 7,, which are known by 
now. We might expect the next approximation to c,(x) on the surface of 
the airfoil not to coincide with the equivalent result obtainable by potential 
theory. 


6. Second approach: flow behind plane shock wave taken as zero- 

order approximation 

For high values of M, and of 5, our method can be improved by assuming 
as zero-order approximation the flow behind a plane oblique shock wave 
generated by a wedge of angle equal to the airfoil tip angle. 

It would be interesting in this case also to try to obtain approximations 
to the flow in the whole disturbed region by expressing airfoil shape, shock 
profile, and perturbation stream function by (23), (24), (25), and (26). But 








AN APPLICATION OF CROCCO’S STREAM FUNCTION 13 


since now the vorticity term in (10) gives a contribution in the first-order 
approximation, the main difficulty of the problem will not be removed by 
the above expansions: the first-order approximation to the differential 
equation will contain a term which depends on the first-order solution itself. 

In this case we shall therefore investigate only the flow in the vicinity 
of the tip by expanding the airfoil and shock-wave profiles in Taylor series. 
Write 


"(0 "0 
y = g(x) = CO) Tae... (52) 
for the airfoil profile, and 
7 = Misin*v = y)+7, ¥+72 '¥?+-... (53) 
or 7+ = tan(v—85) = t9+7,2+722?+.... (54) 


We consider only a region of flow in the vicinity of the leading edge, small 
enough so that it is justifiable to assume 


ny Vi = O((n, ¥)‘) and 1,2* = O((r,2)'). 
The constants 7p, 1, 2) To» 71, and 7, are related by the equations 
No = Mj sin, T) = tan(vy—5Sp) 
m1 = G(M,, 89) (55) 
ne = H(M,, 8973+ L(M, 89)72 
with 
he 29 C08*(vg—8o) 1 ' 
e tan vg tan(vg— 8p) U,(1—uz)y-» 





H = ||_?” _ tan%,—3 tan», tan(vy—8))|>— 
~ 4\ly-1 ’ aa 


(56) 
2y 1 2 1 
(y—1) {no+2/(y—-1)}  (y—-1) ese! 
G tan vg 
2 149 €08*(v9—8g) 








We shall also expand the perturbation stream function 


b = dy +p.t+yH5+... (57) 
and shall assume tentatively 
$; = OY), 
and tie = Om), Py = Om ¥). 
It will be seen a posteriori that all these assumptions are consistent. 
It should be borne in mind that according to our convention the z-axis 


is now chosen in a direction tangent to the airfoil at the tip, and not in 
the free-stream direction. 
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The expression P(w) in (10) can be expanded by (12), (13), and (57) into 
a power series in ¥,, and y,, 





- P(w) = Foot Pro tiyt+ Fo fizt--- (58) 
wil 
1 
sl v(y—1) a pes 
Fy, = 0 
Also Qn) = An oft +o 1¥+f (60) 

- QAno) = Sr 

noino— (y—1)/2v no +2/(y—V} (61) 

O(n) _ 2 1 1 l 








QW)  m—-1 1 M—(y—V/2y m9 +2/(y—}) 
If all the expressions in (10) are similarly expanded by (12), (13), and 


(57), and terms of equal orders are equated, the following differential 
equations are obtained for the first and second approximation: 


Pryy— Bree = Puy Q(9)G7, (62) 
Yoyy— Pore = Py W(9)G7, Pry + 


+ Pro celica ay— (63) 
o) ° . 
BY 


4 on B+ 
—(y-4 1) p2 oe (6? + Debye Prey — —1)— Bp a 


And in general, for the ith approximation, a differential equation is 
obtained of the form 


Piyy—Bbizsr 
— Fy 2, Pry Vizxs Prey» Viyy grees Pi-ryy’ Nor Mas-++> 1:3 Wi-2)- (64) 
When all the approximations up to the (i—1)th are solved, the right- 


hand side of (64) is a known function of 2 and y, and we are faced with 
the problem of solving a differential equation of the form 


Pyy—BWez a d(x, y) (65) 
subject to certain boundary conditions. The general solution of (65) may 
be put into the form 

y 2+By-D 


l ral 
Wx, y) = m(x+By)+n(2—By) +55 | dt | $é,t) dé, (66) 
as can be checked by differentiation. . 
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As for the boundary conditions, the approximate tangency conditions 
are obtained from (16) by expanding around y = 0 and substituting for 
g(x) from (52). For the first approximation the condition 

Py 2(2, 0)+9"(0)x =0 (67) 
is obtained, and for the second approximation 


axl, 0) +9" (O)x pyy(a, 0)+49"(0)a*dr2y(2, 0)-+49"(0)2® = 0. (68) 

The approximate shock-wave conditions are obtained by expanding 

(21) and (22) using (12), (13), (14), and (54), and by evaluating quantities 

on the shock profile by their expansions around the plane y = 7,2. The 
zero approximation shock conditions are 





ae = 1, (69) 
Ug Uy 
: ct tan Vo { y- 1 
(y-+1) uy tan(vy—3)——~ tan v9), 70 
(y+1) ug 1+ tan*(v,—d,)| an(vg—5o) nna an Vp (70) 


the first approximation is given by 


Athy + By py2tQne = 0, (71) 
Agty + Bypi,+O,7,% = 9, (72) 

and the second approximation by 
Ay boy + By tort, 7227+ D, 5 2* = 0, (73) 
Ag oy + Bo por+Cz722*+D, 7} x* = 0, (74) 


where the derivatives of %, and ys, are evaluated at y = 7,2 and the co- 
efficients are given by 


A, = —tan*p, ) 
B, = —tan(v_—8o) 
C, = cos?(vy—8,){tan vp—tan(vg—5Sp)} 
2 
D, = -tan'jg 20 — tan yy—By) 22+ (Pt) + 
27; 27; 12, 
Yaz)” —8,){ Yaz) Yaw) _ (Hae 
+ q(t + by tan(ro—8o) 7,2)\7,2] \r,% 
A, = 0 
B, = 1 
4(v,—8)(3— 
CO, — 008" Yo—%)(8—7 ey hal 
of amie = aay vp+tan(vy—5,) 
(75) 
—tan v, tan(v,—5,)[tan vo tan(vy—3y)]] 
9 2 1 
= ——— cos*(v,—5, pons ~ TE Pe 
(y+1) oN 08%(r9—3y) * Mo 
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D _ [Yazy 4 1 +8 tan vp tan(vy—8,)cos*(vy—8) (4) sh 
a tan vy 7,2) 


- by ( Y2|(S) 4+ ——- = "(% —>o) {1+ tan vp tan(v»—5)} 


tan y, Vo 





27; 


1 
" iv ip pi nro tant »—8,) + tan(ry—85)} — 





9 
— —— cos*(v,—8,) tan 
y+1 (¥9—99) Yo } 


Here again the derivatives of y, are evaluated at y = r)2. 


7. The first approximation 
The general solution of the first approximation differential equation 
(62) is found by application of (66) to be 
Py(x,y) = m,(x+By) +n, (e7—By)+ dP Qn) Gri y*. (78) 
The functions m, and n, are determined by the shock-wave conditions, 
(71) and (72): m,(2) = }.M,7, 22, 
m,(x) = phere 
Pe (A;- —«t,A,)—BA, 
: er . 
as —(As—K79Ay)—Bdg 
2BA,(1—fro) 


with 








where (14) 


cos*(v,—5y) 


x = P,Q(n,)G = 
00 2M tan?y, tan v, tan(vy—55) 





[tan vy — tan(ve—,) |’, 











A, B 
A, =|",' 1) — —tan? 
| 4, B,| 
A, C 2 cos? l 
A, =|~,? 1) s(y,— 2 Me , bt 
‘ A, ©, y+1 ee ral sate Ba | 
aie od 
P B, C,| tan vy 


x tan’) +-tanty—By) —2 2th tan vp tan(v»— >)| : 
y+! 


Hence we may write 


, = (Ax*+2Bry+Cy?)r, (77) 
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with 
A = 4+) 


aR? cos*(v,—5,) tan*(v,—8,) 4! care my 
y+1 {1—tan*(v,—6,)/tan*y,.}| tan®p, eet 
B= $B(4,—) 








2 8in(v¥y—5,)cos*(v_g—5p) 
y+l tan*y9{1—tan*(vp—3,)/tan* mit . Ament a 
C = BUM+M) +46 
Shee Bae) A cos*(vy — 5)[tan vy — tan(ve— mu 
tan*y,| 2 tan v, tan(v)—dq) 








} 

(78) 

The tangency condition (67) gives then the relation between 7, and 

g"(0), which determines the ratio of curvatures of the shock wave and the 

profile of the airfoil at the leading edge. Indeed, the curvature of the 
airfoil at the leading edge is given by 


K,, = 9"(0) 
and the curvature of the shock-wave profile at the leading edge by 
| a" (dr/dx),_» 





sh, = (1-73)! = ©08°(v9—8,)7}. 
Hence, by (67), 
Ken, __ _ ©08*(vg—5p) 
K., 2A 
(y+1) 1—{tan*(v,—5,)/tanyu,} 





= ’ » (79 
4c08(v)—5,) {tan*(vy—8,)/tan®2.} + ${cos*u,/cos*(vyp—3) + 1/n9} (79) 


which is a well-known result (see Kraus (12)). Other equivalent formulae 
have been obtained by Thomas (10), Cabannes (11), and Schaefer (8). 

In a similar way, using the first-order results, (77) and (78), in the first- 
order approximation to c,, the corresponding relation for (dc,/dx),., in 
reference (12) is obtained. 


8. The second approximation 


The second approximation differential equation (63) can now be deter- 
mined by the substitutions 


ty, = 27,(Ax+ By); ty = 27,(Br+Cy) 
ire = 27,4; Prey = 27, B; Pryy = 27,C. 
The resulting equation may be put into the form 
boyy —B ore = % TEX +(%g Ti +03 72), (80) 


5002.41 Cc 
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where a, = 2B{(D—4(6?+1)A} 
a = 2CD—8(6?+1)B°+EF (81) 
LA tan v9 
7 No coor, —8) 
with 
_ 2(y+1) (1—tan*y,)e08*(v)—5y) (—1)? ~ 
Y(y—1) — tanvytan(yy—5 9) = {no—(y—1)/2v Hino + 2/(y—1)} 
—2(1+-tan*y,){(y+1)A+(y—1)C} 
peg \. Oe CO8*}12 CO8*(Y¥y—Sy) no(o— 1)? 
y(y—1) tan*v, tan*(v9—8p) {no>—(y—1)/2y}{ no t+ 2/(y—1)} 
Feat —4tan*v,— tan vo tan(yy—Bo)} + 
ys No 
2 2y—1) l alee l 
mo—) (y—1) fm +2/y—D}_ (Y—D) fm —G— D227} J 
(82) 


The second-order shock-wave conditions (73) and (74) are determined 
similarly by the first-order solution. We obtain for D, and D, 


D, = —C tan*u,— B tan(vg—8,) + 4a99(B+7, C)?+ 
-+-4a9.(A +7, B)?+-4b,, tan(v,—8,)(A +7, B)(B+7,C)— 
—2(A+7, B) 
D, = B—4b,,(A+7, B)(B+7,C)+ 


9 0087(v)— 


8 
tne 0) (A +7, B)[1+3 tan v, tan(vy—59)|+ 


. (83 
c08*(v9—5p) oe 


2 [1-+-tan vp tan(vy—5,)| x 





tan vg 
x Fem aT Vy tan(vg—5o) +tan%(y,—8,)| om 
yt+tl y+l 
9 
aap: oS cos*(vy—5,) tan vo 
y+l 


Finally, the second-order tangency condition (68) becomes, after substi- 
tution of the first approximation terms, 


Yox(x, 0)-+[3 Br, 9"(0)-+49" (0) x? = 0. (84) 
The second-order shock wave and tangency conditions, as well as the 
expanded expression for c,,(x), involve only derivatives of y,. Since we are 


mainly interested in determining c,,(x) to the second order, it is not necessary 
to calculate ys, by (66). It will be sufficient, instead, to determine ¥,, and 
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ey. Using $(x,y) = ay 7Ex+ (aq 7f +05 72)y, (85) 


the following expressions are obtained for these derivatives: 


| e 
Yue = mi(x+ By)+ni(x—Py) + 2 (86) 


Poy = i m3(x+By) —n(x—By)| +o 7 LY + (a Ti+, T)y” 
As in the case of the first-order approximation, so in the second approxi- 


mation, the functions mj(x) and n,(x) are determined by the shock-wave 
conditions (73) and (74): 


























[(oy tot Hotg 78-+Ag/Ay)-+ (Bony 73+ Ag/Ay Br 3 + 
m; (a) = ——— +[( ($a 75— A,/A,)+BA,/A, }r 2 
: 2B(1+ Bro)? vies 
[ (a T9+ ba 7§-+As/A,)—($0 7§+A,/A, \B}ri+ 
n,(x) = +[(4573— A,/A,)— —BA,/ Ail. x2 
. 2B(1—Br»)* j 
where 
_|4, B _|A, G, AR 6 
Ras a a fem le ak She @ po 
=|) D, A, = 
D,|’ . B, D, 











Here 7, is still unknown and must be determined by the tangency con- 
dition (84). The substitution of the above expressions for m}(x) and n,(x) 
in (84) yields the following result for 7,: 











Z,(My, ¥9)72+ Z_(M,, v9)73 = 9” (0) (89) 
with 
, 1{ (dog T5—Ag/A,)+BA/A, — ($a375—As/A,)—BA,/A, 
. B\ (1+ Br»)? (1—Br,)? 
1 (xy To + dorg 78 —As/A,)+( $a, 7§+-BA,/A,) 
4,=3 _ 90 
a (1+ Br) o) 
_ (% To bug T§—As/Ay)— ($e 78+BA,/A,) 12A 
(l—Bro? }+1248 





The derivatives ¥,, and #,, are thus completely determined by (86) to 
(90). 
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A second approximation to the pressure coefficient can now be evaluated. 
By definition 
. (? 1) (91) 
Cc, = -=s| — — . 
> yMi\p, 


But since the flow along the airfoil is isentropic, 


pare Pa( ‘iy = gil 


Py PNT ~ p,\l—u3 





(92) 


This expression can be expanded by (12), (13), and (57), and the values 
of the derivatives of %, and ¥, on the airfoil surface, which appear in this 
expansion can be further expanded around y = 0. We obtain, to the second 
order of approximation, 


, —1) 2 
=)" 1 = $4 YM ya, 0)+ 


1 —u3 Mi— l 


M3 —1)M$+ M2+1 
+ xa ae 1 ee Piy(z, 0)+-y7,(x, 0) 4- 





4 aba (2, 0)-+29(Wayy(2, 0). (93) 


The functions %, and , determine the first and second derivative of c,,(z) 
at the leading edge. 

A check of the results can be obtained by calculating the approximate 
shape of the shock-wave profile by 


T T : 
y = mr+ 5a*+ got 


and comparing the result with the curve obtained by characteristics (4, 15). 
Fig. 1 shows this comparison for a 10 per cent thick parabolic biconvex 
airfoil placed symmetrically in a uniform stream at M, = 00. In the figure 
are plotted the shock-wave profile obtained by the method of characteristics, 
the results of the present first and second approximation, and the curve 
obtained by the shock expansion method (4). The second-order curve is 
closer to the curve obtained by characteristics even relatively far from the 
leading edge. 

Fig. 2 is a plot of the results for c,, for the same airfoil and flow. In Fig. 3 
are plotted values of c,, for the same airfoil in a uniform stream at M, = 10. 
Our second approximation curve coincides with the curve obtained by 
characteristics (15) up to z/c = 0-35 within an order of 2 per cent. 
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SUMMARY 


A two-dimensional jet of incompressible, inviscid fluid is deflected by an aerofoil 
immersed in it. The relation between the aerofoil incidence and jet deflexion angle, 
the extent to which the aerofoil displaces the jet from the position it would occupy 
in the aerofoil’s absence, and the ratio of the fluid flowing above the aerofoil to that 
flowing below are calculated in this paper. Other results, such as the velocity 


distribution over the aerofoil and the shape of the free surface of the jet, are also 
given. 


1. Introduction 


FREE jets of air and water are frequently used to test the performance 
of aero- and hydrofoils, and an ‘exact’ theory for dealing with this type of 
flow is of some interest. 

One application of such a theory is to flexible walled wind tunnels (1). 
In these tunnels the walls are first adjusted in shape until the pressure on 
them is constant, and the flow then corresponds to that of a free jet about 
an aerofoil (neglecting the effect of the boundary layer on the walls). The 
walls are next moved according to a rule—which must be derived from 
theoretical considerations—until they lie along the streamlines appropriate 
to an infinite stream. The existing theory for a lifting aerofoil is based on 
the theory of a vortex in a free jet (1); this could be improved with the aid 
of the theory given in the present paper. 

Havelock (2) has given a theory which applies to a thin aerofoil between 
free surfaces, but in this theory the free surface boundary condition is 
applied on lines parallel to the y-axis instead of on the free surface itself. 
But even the assumption that the aerofoil’s incidence—and hence the jet 
deflexion angle—is very small, must still mean that downstream of the 
aerofoil the jet diverges further and further from the y-axis by an amount 
which increases infinitely. There is some a priori doubt, therefore, as to 
whether Havelock’s solution is relevant to the free jet and aerofoil problem; 
in any case the validity of his approximate boundary condition ought to be 
checked by a theory in which the existence of a finite deflexion angle is 
taken properly into account. 


+ Nuffield Research Professor in mechanical engineering. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958) 
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It is found in this paper that (see section 11) even in the limit as the jet 
deflexion angle tends to zero the author’s results are not in agreement with 
those given by Havelock. This is found to be due to a discrepancy between 
the number assumed by Havelock to be the incidence and the real incidence, 
and a simple modification to Havelock’s incidence brings the results into 
agreement. 

The main difficulty of an exact theory of the aerofoil and jet problem is 
that the position of the free surfaces is not known initially, and furthermore 
this position is dependent on the circulation about the aerofoil. It is there- 
fore necessary to assign a value to the circulation from the beginning, but as 
the theory is developed sufficient relations emerge to enable the circulation 
and other similar parameters to be eliminated in favour of geometrical 
variables such as the incidence, jet width, and chord length. These relations 
are summarized in section 10, and readers concerned only with the applica- 
tion of the theory can turn immediately to this section. 


2. The conformal transformations 


The transformations and planes used in the subsequent theory are shown 
in Fig. 1. The z-plane (z = x+-iy) is the ‘real’ or physical plane; in this 
plane the aerofoil (A BC B’A’) is shown deflecting a free jet (EF, GF’, E’) 
through some angle 8, say. There will be a circulation about the aerofoil, 
so in the w-plane (w = ¢+-%) of equipotentials (¢ = constant) and stream- 
lines (ys = constant) it is necessary to have a jump in potential across some 
barrier such as AZ or B’ F’,, which reduces the connectivity of the region. 
AE is the most suitable choice in the general case, for then the w-plane can 
be arranged symmetrically as shown in the figure, and this in turn simplifies 
the subsequent transformations. The origin in the w-plane is placed at C, 
a point midway (in the w-plane) between the front and rear stagnation 
points B and B’ at ¢ = —2a, 2a respectively. If one side of the barrier, 
A'E’ say, is at ¢ = 41, and the other side, AE, is at 6 = —4I, then the 
circulation will equal I. 

From the point of view of potential theory the shape of the w-plane 
shown in Fig. | israther awkward. A more convenient shape is obtained by 
mapping the w-plane into the rectangle shown in the {-plane ({ = y+in). 
The t-plane shown in the figure is simply an intermediate mapping used in 
deriving the relation between the w- and {-planes. As each of these planes 
are (generalized) polygons in shape, each can be mapped on to the t-plane 
by the Schwarz—Christoffel mapping, and then eliminating ¢ from the two 
resulting mapping equations yields the relation between the {- and w-planes. 
The following table gives the values of w and { at corresponding points, and 
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together with Fig. 1 serves to define all the constants introduced in the 


mapping. 
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In terms of the symbols defined in this table the relation between { and w 
is readily found to be 


w = * (hy +hg)[{2Z(6)-+ZA—8)—Z28+.2)}L— 211, 3)], (1) 


where Z(x) is Jacobi’s zeta function, and II(z,a) is Jacobi’s incomplete 
elliptic integral of the third kind (see (3), pp. 518-23), i.e. 





aa ( k*snacnadnasn*x 
aa tad | 1—k*sn2a sn2z - (2) 


0 
The real and imaginary quarter-periods of the Jacobian elliptic function 
involved here are the numbers K and K’ introduced in the table above, 
while & is the modulus of these functions. 

We shall assume for the moment that a, [’, h,, and A, are known quantities 
(it will be shown in section 7 how these are determined) so that our mapping 
into the (-plane introduces three unknown constants, viz. 5, A, and k 
(K’, K are determined uniquely by k). Relations defining these constants 
are found by substituting the corresponding values of w and ¢ at points 
A’, B’, and E’ in igs (1). It is found that 


ik = {2Z(8)+ Z(A—8)— Z(6+A)}A—211(A, 8), (3) 


aT 
2(hy +hg) 


and 23 = eae 


= {Z(A—8)—Z(8+A)}K, (4) 


'K’ 


K 5 
ithe 7} Ti ak (9) 
In the next section we shall also need the result 


dw _ 2th hy) sndensdn8({ sn?A—sn®Z ¢ 
a ri?" 1—sn% sn?A | 1 —k*sn%3 ai 


which follows from (1), (2), and the relation 


2k? sn*A sn8cn8dn5d 
1—k*sn?8 sn?A 
The importance of the {-plane is that it arranges the mixed boundary 
conditions of the flow problem in a very simple pattern—the known aerofoil 
shape lies on » = 0, —K < y < K, the constant pressure surface is on 
» = K', —K <y < K, while any function f(y, 7) characteristic of the flow 
satisfies f(—K,») = f(K,n), 0 << K’. We shall therefore adopt the 
{-plane as the plane of independent variables. 
Before introducing the dependent variables we mention here one impor- 
tant special case of (1) that has been used previously (for example, (4) 








(6) 


2Z(8)+Z(A—8)—Z(8-+A) = 
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and (5)). This is the case in which the circulation is zero, i.e. ! = 0, and the 
stagnation streamline y = 0 lies midway between the free boundaries of the 
jet, ie. hy = hy. With these values (5) reduces to 8 = 4K, then (4) shows 
that A also equals 3K. Putting a = 4K in (2) and making use of Landen’s 
transformation (see (3), p. 508), we find that 
. 1+k,sn{(2K,/K)l,k ) 
I({,4K) = —Z(4K)C—4fln A A ah 7 

(9%) Ges (i —k,sn{(2K,/K)¢, ky} "7 

where the modulus of the elliptic functions appearing here is 
1—k’ 
= 1+k” (8) 

k’ being the complementary modulus (1—k*)!. The corresponding quarter 
period K, is related to K by K, = 4K(1+4’). Equation (7) permits (1) to be 
rewritten 





tanh =k, anh) (9) 


where h = }4(h,+A,). Similarly (3) reduces to 
ky = tanh. (10) 


The importance of this special case lies in the fact that the stagnation 
streamline B’ F’, now lies on the straight line y = 4K, so the transformation 
{’ = (+4K shifts B’F’ to (’ = +K. In the general case the stagnation 
streamline cannot be placed on opposite sides of the rectangle in the ¢- 
plane in this way. Having the trailing edge streamline on (’ = —K-+0, 
{’ = K—0, enables us to deal with small unsteady perturbations of the flow, 
for such unsteadiness produces a vortex sheet of calculable strength lying 
along the stagnation streamline, i.e. a jump in velocity from ¢’ = K—0 
to ¢’ = —K-+-0. And this type of boundary condition can be treated (5) 
by an extension of the theory for steady flow given in the next section. 


3. The general theory of the flow of a jet past an aerofoil 
The only satisfactory dependent variable to use in the theory is that 


defined by y . 
f= In(e) = n(Z} +40, (11) 
dw q 


where (q,@) is the velocity vector in polar coordinates, and U is the (con- 
stant) velocity at the free surface of the jet. The angle 6 will be measured 
from the positive direction of the z-axis. By (1) and (11) the complex 
function f is an analytic function of z, w (as w = w(z)) and of {. Thus, 
in the rectangle in the (-plane we have to find an analytic function f(Z) such 
that imf = 0,, the known slope of the aerofoil surface, on 7 = 0, ref = 0 
on the free surface » = K’, and {(—K,7) = f(K, 7). 
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The solution to the mixed boundary value problem just defined is known 
to be (6) 


K 
" (K, 
100) = | oS (RO —Ds ka] dy* 
-K 


=} f 0, d{in| sn(32(y*—,,)]], (12) 
ee 


where K,, K, and k, are the numbers defined in section 2, and by sn’/sn(w, k,) 
is meant (d/du)ln{sn(u, k,)}. (The result given in the author’s earlier paper 
is actually the special case of (12) in which K = 2K;,.) 

A more convenient form of (12) can be obtained by applying Landen’s 
transformation. As K,/K = }4(1+k’), it is found that the transformation 
enables (12) to be written 


K 
Ao) = = [On dfln(snfi(y*—2), Hed{H(o*—D), #)} 
=i 


l 


= se J Padtin(an'{y(y*— Cea Hey*—D)} 


2n 





K 
l fin 1—en(y*—{) 1+ en(y*—2) 
= — 6,d 
ae | aie terns *—L) 1+dn(y* al} 


where the modulus of the elliptic functions k is omitted for brevity. The 
modulus of all elliptic functions appearing below is also k, so there is no 
ambiguity involved in the abbreviated notation sn u, etc., which will be 
used in the remainder of this paper. 


Performing the differentiation in the last expression for f({) we arrive at 
the concise result 


K 
HO) =~ | Agcsty*—2) dy*, (13) 
—-K 


which is the required general solution of our boundary value problem. 
4. The jet deflexion angle 

Let the jet be deflected through an angle f in the sense shown in Fig. 2, 
and let the z-axis be selected so that upstream at infinity the flow is parallel 
to it, ie. Os... = 0. Then _,.. = —B. 

On the free surface { = y+7iK’, and on substituting this in (13) it is found 


that the slope of the free surface is given by 
K 


Ay) =~ | O,dn(y*—y) dy’ (14) 


—K 
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The point upstream at infinity is at y = —é5 (cf. Fig. 1), and at this point 
@ = 0, so that from (14) we have 
K 


ut | 9, dn(y* +8) dy*. (15) 
” ie 
Similarly, 6(5) = —B, so that 
K 
p=—> | 0, dn(y*—8) dy*. (16) 
7 
-K 


Equation (15) serves to define the position of the front stagnation point 
relative to some fixed point on the aerofoil (e.g. the leading edge), while 
(16) determines the jet deflexion angle. 
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Fic. 2. The jet displacement 


5. The closure conditions 


No restriction has so far been placed on the magnitude of the deflexion 
angle 8. However, a restriction is necessary if the relation z = z(f) is to be 
expressed in a tractable form. Then the closure conditions, i.e. that 
$¢ dz = 0 about any circuit enclosing the aerofoil, can also be expressed 
simply. Thus, for mathematical reasons we now restrict 8 to be of first 
order in magnitude, and ignore second and higher order terms. From a 
practical point of view this is not an important restriction, for the aerofoil 
incidence « must always exceed the deflexion angle (see equation (44) below), 
and the usual restriction on a, i.e. that it be less than the stalling angle, will 
also restrict £. 

If the z-plane origin is at C (Fig. 1) where w = { = 0, then (11) and (6) 
= 








dt. 
k?sn8 sn?A 1—k*sn*5 sn? $ 
. (17) 


F 2) en2 
Us = | esos epee = 2 hth) ahattel f y _sn?A—sn®f 
0 
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The closure condition, namely $ dz = 0, can be written (see Fig. 1) 


K+in d 
w 
| | so (Tr) dt =0, (18) 
—-K+in 


where 7 is any number between 0 and K’. On the jet surface and at large 
distances from the aerofoil it follows from (11) that f = i@, where 0 < B, 
so if a value of » nearly equal to K’ is used in (18), we can write ef = 1+-/f to 
first order in 8. Then as w(f = K+in)—w(f = —K-+in) equals the 
circulation I’, (18) becomes 





K+in d 
r+ [ s0(q) a =o. 
—K+in 
or by (13) and (6), 
K 
P42 | O(1K-+in,y*)—(—K-+in, y*)} dy’, (19) 
where SK t 
2 k*sndenddné [ cs(y*—{)(sn?A—sn*f 
Hy") = =th+h) a peienh | T= Beer 20) 


0 
This integral can be evaluated by using the formula 
_ snfen{dny*+sny*eny* dnl 


ces(f—y*) sn*y*—sn®{ 


It is found that 





k*sncnddné 
1—ksn% sn®A * 
{ sn*y*—sn?A 1—dn(y*—f) 1+dn y* 
. \1—k*sn*5 sn*y* =| 1+cdn(y*+) 1—dny* 
2(1—k*sn%$ sn?A)/sny*eny*, ok 
“Kat l aso 
__dny* 
k*sn dnd 
On taking account of the periodic properties of the various functions in 
this expression, we obtain 
I(K+in, y*)—1(—K-+%n, y*) 
k*sn$ cnd dnd - 
(1 —k*sn6 sn?A)(1 — k*sn*8 sn*y*) 


Ig, y*) = = (hy +h) 








{tan-Y(dn {808)—tan-Y(s08)}]}. (21) 





1 
= = (hy +h,) 
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sn y*cn y* 


x | 2mi(sn*y* —sn?A)+ 2n(1 —Kan% sn?A) dnd 
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Substituting this into (19) the imaginary part gives 








K 
sn*y*—sn?A aE. 
| 9 1—k*sn*5 sn*y* dilaton = 
and the real part gives 
1 2k*sn § cn sn y*en y* 
ed dy* = 0, 
P+ oh +hy) { 9 i—sn%sn4y* ” 
-K 
, K 
or P—2 (hy +h) | 9,{dn(y-++8)—dn(y—8)} dy = 0. 
-K 
Applying (15) and (16) to this result we get the concise relation 
TP = (h,+h,)B = HUB, (23) 


where H is the width of the jet at infinity and U is the corresponding 
velocity (h, +h, is the total change in the stream function across the jet, and 
therefore equals HU). 

Equations (22) and (23) are the required closure conditions, but it must 
be remembered that they are valid only for small £. 


6. The relation between z and ¢ 


Equations (13) and (17) together give the exact relation between z and (, 
but, except for some special cases, the integration must be performed 
numerically. The approximation ef ~ 1-+/f can be applied to (17) only if the 
incidence is small, and if the path of integration does not terminate near 
the leading or trailing edges where |f | is apt to be large. This restriction 
on « is more stringent than that imposed in the previous section, as f is 
always less than a, and for a wide jet is much less than «. 

When the approximation is valid it follows from equations (13), (17), 
and (20) that : 

Uz =w+- [ 6, 1(L, y*) dy*. 
"ek 
Now equations (15) and (16) may be combined and written in the form 
nm _2 fsacee POT... [ Phreratete 
7 1 —k*sn*6 sn*y* 
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so that when (21) is substituted into the equation for Uz the result can be 
written 


Uz = w— ' (hy +hy)B{tan-(se { dn 8)—tan-(dn {se 8) -+-tan-1(se8)}+ 


K 
. —_— 
+g 4a) Seae ’ sn*y*—sn?A 
td 1—k*sn*6 sn?A 1—k*sn*6 sn*y* 
-K 
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7. Position of the aerofoil within the jet 

The actual position of the aerofoil relative to the jet will obviously depend 
on both the aerofoil incidence a, and the position of the aerofoil centre 
relative to the jet position upstream at infinity. Let JK in Fig. 2 be the 
centre line of the jet in its undisturbed upstream position, and let the 
vertical distance from C (origin of the z-plane) to JK be D. Then if D is 
fixed h,/h, will clearly decrease with increasing «. The relation between 
these variables can be found from equation (24). 


Let « be a small positive number, then lim im 2({ = «—5+7K’) is the 
«e—0 


value of y on the upper free surface upstream at infinity. The jet centre- 
line JK is a distance $H below this free surface, and it therefore follows that 


D = lim im 2(f = e—8+1iK’)—4H, (25) 
«0 
where z is given by equation (24). We shall take the limits of the terms in 
(24) separately. 

First, combining equations (1), (+), and (23) we find that 
d,{(7/2K)(C+8) |(¢K'/K)} 
54{(/2K)(l—8) |(¢K"/K)} 
where 2{II(¢,5)—{Z(8)} has. been replaced by the logarithmic term 
involving the theta functions (see (3), pp. 479 and 523). Thus, from the 
quasi-periodic properties of the theta functions, 


lim w(f = «—8+iK’) 
«0 
_ Blhy+h,) ((K’—8) | 1 —e8iK 9 (me/2K) 
- 2K +; 5,(n3/K) |’ 
the parameter of the theta functions (r = iK’/K) being omitted for brevity. 
Equations (5) and (23) permit this limit to be expressed as 
8, (me/2K) 
5,(73/K) 


(26) 
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lim w(t = ¢—84-iK") = thy + © (y+) im | (27) 
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Secondly, 
lim {tan-'se(e—5+¢K’)dn§—tan-'dn(e—8+7K’)sc 5+ tan-"(se 8) 
«0 
i, k*sn*6 en 
‘ie 5in a tan-!(csd), (28) 


and, finally, 


ot 1—dn(y*—{) 1+-dny* sae 1 x 9 
lim im In( = +dn(y*—t) 1—dny* ) = 2 tan~les(y +8). (29) 


From (23) to (29 





9) the following expression for D can now be derived : 
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2D _ h,— B,., k*sn*6.cn?3] _ 
H =e 7 dn*5 


K 
__ 4 Ksndcnddnd sn*y* —sn?A 


nm 1—k*sn% sn?A “1—k*sn%5 sn*y* 
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tan~'es(y*+-5) dy*. (30) 


Another important geometric number is the chord to jet-width ratio, 
c/H. A simple expression for this can also be derived from (24) provided 
it is assumed that the distance between the front and rear stagnation 
points is a sufficiently close measure of the chord length (always true for 
aerofoil shapes at moderate incidences). In this case 


¢ = 2(f = A)—2(f = —A) 
and so (24) yields 


c “F tan- basil = 1 Hen anfh 


H ~ h+h,* 
K ‘ 

x fo. Smty*—sn*A_},,|1—dn(y*—A) 1+dn(y*+A) 

| “| —#sn% sn®y* | |1-+-dn(y*—A) i—dn(y*+a) 





dy*. (31) 
-K 

It will be recalled that in section 2 we assumed temporarily that a, I, 
h,, hg were known quantities, but of course in a practical problem of an 
aerofoil in a jet the known quantities are most likely to be c, H, D, and «a. 
However, we now have sufficient relations, viz., equations (16), (23), (30), 
and (31), to calculate a, T’, h,, and A, from given values of c, H, D, and «. 
(The incidence « is implicit in the right-hand side of (16), see section 9.) 


8. The lift and circulation 


The equation of linear momentum shows that the force on the aerofoil .L 
makes an angle of 4(7—f) with the y-axis, and that its magnitude is 
2HpU* sin $8, where p is the fluid density. Thus, to first order in 8 we have 
L = pHU*8, and it therefore follows from (23) that L = pUT, in agreement 
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with the classical result for an infinite stream. The lift coefficient, 


C, = L/(tpcU?), 
is given by C, = 2HB/c. (32) 


9. Thin aerofoil theory 

With thin, ‘unflapped, aerofoils we can assume that 0, = —a, except 
near the leading edge. At the leading edge the flow is réversed in direction 
(i.e. 6, = m—«a) between the stagnation point at = —A and the nose at 
{ = —A-+e, say. Thus, we can write 


= —« (-K< ¥<-—A, —A+e< YS K) \ (33) 


m—a (—A<y<—A+¢) 


where at moderate incidences « will be a small quantity of the same order 
as a or f. 
From equations (15), (16), and (33) we obtain 


B = —e{dn(A+5)—dn(A—8)}) 


a = edn(A—8) j (34) 
to first order in «. Thus 
p= oft pee (35) 
Similarly equations (11), (13), (33), and (34) give 
= = U{1+ia+acs(d+f)nd(A—8)}, (36) 
and in particular, the velocity on the aerofoil surface is given by 
h = 1+acs(A+y)nd(A—8). (37) 


Equation (22) is satisfied by (33) correctly to first order in «. 

In the evaluation of the integral in (24) it will be found that the leading 
edge term in (33) contributes no first-order term. The integral arising 
from the incidence term is most easily integrated by first differentiating it 
with respect to ¢, and then reversing the order of the two integrations now 
necessary. The final result is 





Us = w(1—ia)4- at MMB?) ftan-1(dn { 50 8)—tan-1(s08)}— 


2 Crt MP tan-1(s0f dnd). (38) 
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The equation corresponding to (30) is most directly obtained from (38) 
by making use of (27) and (28). It is found that 

2D h,—h, By, | Ldn 26 4 20), 3, (735/K )dsd ned 

H  h+h,' am |1+dn23| ° a (8kk’K |x)! 


k?sn8 en25 1—dn 28 








(39) 


as 





ans = 1pdn 38’ 
and 
lim 9,(5) /« = 7 84(0) = 57-94(0)0,(0)9,(0) = (2k Km) 


Equation (39) can also be written in the more convenient form 


2D ih, B,/t—tes 25)9,(n5/K)] | a a aS 


H~ hth,’ |  (8k’K/z)! (8k K/n)i 





(40) 
The equation corresponding to (31) from (38) is found to be 


c 4a 2 
= _—.- — —Btan-(seAdn83), 
hak <> ae capac, 
or retaining only the highest order terms, c/H = 4a/(h,+h,), which is 


equivalent to Po: a (41) 


A similar calculation applied to (26) gives 
4a 2H = in( 5 a{(7r/2 ee 


7 = Shae dent Ace! 


B,{(2/2K)(A—8)} 














(42) 
wm, oe 3 4{(7/2K)(A+8)} 
so by (41) ao a aerate: 2] 
10. Summary of equations : examples 
The principal equations of the theory are 
Rie ae (Z0— 5) —Z(3-+A)}, (43) 
B = a{1—dn(A+8)nd(A—8)}, (44) 
Ae ig tus (45) 
hi, +h, C wt dy 
C2) |dyi(w/2K)(A+8)} 
oe - ' (46 
Hw |8{(7/2K\A—8)} 
and 
2D h,— > + Ein inf = 25)3,(73/K) 20-8), (1+dn ee | 
H  ih,t+h, 29,(0) Ae” 29,(0) , 


(47) 














ON THE DEFLEXION OF JETS BY AEROFOILS 37 


The first of these five equations follows from (4) and (23), the second is the 
same as (35), the third follows from (5) and (23), the fourth and last are 
respectively (42) and (40). 

Suppose that we have a given geometric arrangement of jet and aerofoil, 
then c/H, D/H, and « are known quantities. Equations (43)—(47) fix the 
values of the five unknowns f, K, A, 5, and h,/h,. 

An easier procedure is to assign a set of values to A, 5, and k. Then 8 
follows immediately from (43), « from (44), h,/h, from (45), and finally 
(46) and (47) yield c/H and D/H. After a few trial values of A, 5, and k, it is 
not hard to find values that will very closely produce any desired aerofoil- 
jet combination. By way of illustration we consider two examples. 


Example 1. Aerofoil in the centre of the jet 
If the jet is undeflected it fellows from (43) and (45) that8 = 4K,A = 4K. 


Then (46) gives ce 2, /8,(0)\ 1, /1 
; = ee 3 poe eeepc nits 
a= = In( sia} ms ~In(7)) (48) 


If the aerofoil is now at a small incidence, 5 and A will differ from 4K only 
by terms of order «. Thus, to first order in a, (43) gives 





8 = o(1—dn K nd0) = a(1—k’) = (1—e-*-/#), (49) 
and if c/H is small, we find from (32) and (49) that 
aC, 2H ce |, we wife )' 
a = —e-me SS 4 L.. — —{— . 
Fo aire ies oo 


As dn(A+-5)nd(A—8) is a minimum at A = 4K, 5 = 4K, it follows from 
(43) that (eC,,/@«),_, is a maximum when the aerofoil is in the centre of the 
jet. 

Example 2. k* = 0-9, 8 = 1-00, A = 1:35 

This set of values and equations (43)-(47) yield in turn 8 = 0-1857 
(= 10-64°), « = 02679 (= 15°34°), 

(hy —h,)/(h, +h.) = 0-340, c/H = 0-372, and 2D/H = 0-268. (51) 

Clearly only a small computing effort would be required to establish 


sufficient entries in a table, which would allow interpolation to any desired 
values of a, c/H, and 2D/H. 


11. Comparison with Havelock’s theory 

In order to test the theory given by Havelock (2) it is sufficient to com- 
pare his result for (@C,,/@«),., with that given by equation (50). As Have- 
lock’s results are all given in the form of infinite series (in c/H and D/H) 
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it is not easy to make an exact comparison; he gives an equation which is 
equivalent to 


ac, _ef, wie), # fey c \ 
(=z), = 2['~ Tala) +ip0(@) + (a) } 


and therefore is not in agreement with the author’s theory. Yet there are no 
errors apparent in the theory which yields (52). 

An equation like (52) can be deduced from the author’s theory as follows. 
Replace « in (49) by «’+-48, then by (49) 


1—k’ 7 
= 9n/f |) — Ba’ t athin 
shai Gs) F tanh (Fi) 


and therefore from (32) 


eC, “a oe) .f, «fer, & ic) c \*) 
(Fei) e (Se) = 2! Fale) + ia0(H) + UH) f 


This makes it clear that Havelock’s incidence is not measured from the 
upstream jet direction (as implied in his paper), but from the average of the 
upstream and downstream jet directions. This seems to be a reasonable 
conclusion, for the errors introduced by the misplacement of the free 
boundary condition downstream of the aerofoil must clearly be balanced 
by a like displacement upstream of the aerofoil. 

As Havelock gives no expression for the jet deflexion his results cannot 
easily be modified to allow for the discrepancy in the incidence. 


(52) 
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SUMMARY 


The motion of a sphere along the axis of an unbounded fluid rotating with uniform 
angular velocity is discussed on the assumption that the fluid is undisturbed far 
upstream. It is shown that no wave-like components can then occur upstream of 
the body. The stream function % and the drag are found for various values of the 
Rossby number ka, defined in (1.1) below. When ka = 5-76 both these quantities 
become infinite everywhere and it is inferred that the fluid cannot be undisturbed 
far upstream when ka > 5-76. It is probable that cylindrical flow which must then 
be present must also occur at smaller values of ka of which the lower bound is not 
known. 


1. Introduction 

Tue chief hydrodynamic interest in the motioa of a sphere along the axis 
of a uniformly rotating fluid stems from the unusual nature of the flow. 
The first experiments were made by Taylor (1, 2) who found that at small 
values of the Rossby number 


20a/W = ka, (1.1) 


where ( is the angular velocity of the fluid, a the radius of the sphere, and W 
is its velocity, the disturbance in the fluid was similar to that when 2 = 0, 
dying out as the distance from the sphere tended to infinity. At ka = 2m 
the character of the disturbance became markedly different. Ahead of the 
sphere there appeared a régime, near the axis of rotation and bounded 
approximately by a circular cylinder C circumscribing the sphere an d having 
its generators parallel to the axis of rotation, where the flow was cyli ndrical. 
Inside it the fluid was pushed along in front of the sphere and at the same 
time it behaved with respect to the fluid outside C as if it were solid. 
Somewhat similar results were found by Long (3) using a body with a 
spherical nose and a conical tail. He also examined the flow beh ind this 
body and found that for all values of the Rossby number the fluid n ear the 
boundary passed round to the rear of the body and was not carried along 
with it. 

A theoretical solution when ka is large (4) confirms the general fe atures 
of Taylor’s observations (1, 2) but not Long’s (3) on the flow to the rear of 
the body. The reason for this failure is not known. 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958] 
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Taylor (1) also considered the flow from a theoretical standpoint. Let 
(r, 6, z) be cylindrical polar coordinates of a point P with origin at the centre 
of the sphere and the axis of rotation as Oz, and let (u,v, w) be the corre- 
sponding components of velocity. Further we suppose that the sphere is 
at rest and that in addition to the rotation the fluid streams past it with 
uniform velocity w = —W when z = +00. From the equation of con- 
tinuity a stream function ‘Y may be introduced such that 

Ley 
—-— = &. 
r oz 
Then, provided that the problem is properly posed so that the conditions 
at infinity can be as specified, we have (3) 


op le Ys = 6, 


er? c oz? 


(1.2) 


where y= F+4+4Wr? and rv = —kY, (1.3) 


with boundary conditions 


y>Oasz>oo and & = }Wr’ on the sphere. (1.4) 


Taylor obtained a solution of this equation which included an arbitrary 
constant. Long (3) showed that nearly every solution of (1.2) satisfied the 
boundary condition as z—-00, so that another boundary condition was 
needed to complete the solution. The indeterminacy is due to the exist- 
ence of wave-like solutions of (1.2) which, if the fluid is unbounded, die out 
as r?+-z? + oo. If, however, the fluid is confined within a cylinder of fixed 
radius b (> a) the situation is different. When kb < 3-83 all acceptable 
solutions of (1.2) die out exponentially as z > 00, but if kb > 3-83 there exist 
wave-like solutions which cannot die out as z -> 00. Hence if the problem 
is properly posed so that 4 -—> 0 as z > +00, they must be excluded and 
there is then a unique solution. Long also carried out experiments which 
confirmed that there were no wave-like motions ahead of the body, the 
streamlines only oscillating to its rear. 

Starting from a study of sources in a rotating fluid enclosed within such 
a cylinder Fraenkel (5) obtained a solution to the flow of an unbounded 
fluid past the sphere when ka is small by letting 6 > 00. In adopting this 
limiting process it is assumed that since no wave-like motions occur ahead 
of the body when 6 is finite, neither will they occur when b = oo but, as 
pointed out above, this is not self-evident. The difficulty arises because 
two limiting processes, b > oo and z— 00, are involved and the order in 
which they are taken is important. 
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For an unbounded fluid the most general solution of (1.2), free from 
singularities in z > 0 is 


¢ =r [ A(o)\J,(or)exp{—z(o*—k*)} do (1.5) 

0 
where A(c) is an arbitrary function of ¢. The wave-like part of this solution 
when z is large comes from that part of the path of integration in which 
ao < k where the exponential oscillates. If no oscillations occur when z 
is large A(jo)=0 if o<k. (1.6) 

In section 2 of this paper we show that no matter how the motion is 
started, if ultimately it is nearly uniform a long way ahead of the body, 
y there has the form 

x k 
r { J,(or)A(o)exp{—z(o?—k*)!} do+-r f B(o)J,(or) do, (1.7) 
k 0 

where A and B have to be found. If ¢ > 0asz—> +0, B = Oand ¢ then 
satisfies both the differential equation (1.2) and the extra boundary con- 
dition (1.6). It is quite possible, however, for B to be non-zero, in which 
case the motion is cylindrical at large distances upstream and the problem 
stated in (1.2) and (1.4) is not properly posed. There is thus some theoretical 
doubt as to whether the motion of an unbounded fluid is uniform a long 
way ahead of the sphere for any ka. The experiments are against such a 
cylindrical component to y (i.e. independent of z) for small values of ka, 
Taylor (2), for example, finding that it first appeared when ka ~ 27. It 
must be remembered, however, that they were all carried out with fluid 
inside a tube of finite radius 6 with a/b small, and it has been argued that 
a drastic change may take place in the character of the flow as a/b > 0. A 
reason for not believing this is given below. 

The nature of the flow due to a sphere is examined in sections 3 and 4 of 
this paper assuming that % > 0 as z > +o and (1.6) holds. The complete 
family of solutions y,, (m = 1, 2,...) of (1.2) satisfying (1.6), of which the 
first was found by Fraenkel (5), is obtained in section 3. On writing 


b= LAntm (1.8) 
the solution which satisfies (1.4) in addition may be obtained, the A,, being 
constants which satisfy an infinite set of linear equations. They are found 
in section 4 and are displayed in Table 1 there, together with the correspond- 
ing drag coefficient C, as functions of ka. So long as ka < 2 all the A,, 
except for A, are small, but for larger values of ka they rapidly increase in 
value and the streamlines become greatly distorted. Eventually, at 
ka = 5-76 all the A,, become infinite, because the determinant of the 





42 K. STEWARTSON 


coefficients of the linear equations which they satisfy vanishes. This means 
that the possible values of y on the sphere are restricted and to satisfy (1.4) 
additional terms, which can only be cylindrical from section 2, must be 
added to %. If ka > 5-76 the solution given by (1.8) although mathe- 
matically consistent is unrealistic. From this and a comparison of Cp with 
a known solution, it is concluded that a cylindrical component certainly 
occurs in % if ka > 5-76, probably if ka > 3, but may occur for all ka > 0 
although there is then an apparent discrepancy with observation. It is not 
believed that such a discrepancy would be resolved by considering the 
fluid to be enclosed in a large but finite container for the following reason. 

When ka is infinite the flow, both when a/b + 0 and when a/b = 0, is 
known to be cylindrical. Now let us suppose that when ka is large (but 
< kb of course) ¢ > 0 as z+ +00. Then if kb = 0, from (1.5) and (1.6), 


bw rd(ker) [ Clrje dr (1.9) 
0 
near the sphere, where C is an arbitrary function. On the other hand, if 
kb is large, from the corresponding series for yf (cf. 3), 


gb ~ rJ,(arjexp{—z(a?—k?*)#} 


near the sphere, where ab is the smallest of those roots of J,(x) = 0 which 
are greater than kb. In either case y cannot be « r? on the sphere and 
so a cylindrical term must be present. It is inferred that for sufficiently 
large kb there is a value f(kb) < kb of ka such that if 
ka > f(kb) 

the hypothesis that y > 0 as z > +00 must break down. The way in which 
the breakdown occurs when kb = oo is then likely to be typical, i.e. ¢ is 
equal to the sum of a series of functions corresponding to those in (1.8) 
and at some finite value of ka their coefficients would all be infinite. The 
argument used when kb = oo also applies so that there is a theoretical 
possibility that % has a cylindrical component for all possible ka even when 
kb is finite. The fact that experiment does not bear this out enables us for 


the present to regard it only as a remote possibility both for kb finite and 
kb infinite. 


2. Ultimate conditions far upstream of the body 
In this section we show that if the flow far upstream of the body is 
uniform at all times then (1.6) holds, i.e. no wave-like motions occur ahead 
of the body, no matter how it was started from a state of relative rest. 
Let the body be at rest at a large distance downstream of the origin and, in 
addition, let the motion at an infinite distance upstream (z > +00) consist 
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of a uniform rotation and an axial component of velocity —W. Then, 
however the disturbed motion from this state is started the disturbed com- 
ponents of velocity are already small at z = 0. To describe the motion 
when z > 0 it is therefore legitimate to linearize the equations of motion as 
follows. Using fixed rectangular cartesian coordinates (x, y,z) the velocity 
components may be written 


(—Qy+4,, Or+4,, a W+4,), (2.1) 
where ¢,, 7), 7, are assumed to be small if z > 0. Then if p is the pressure, 
p the density, and P = p/p—4Q%(x?-+y2), (2.2) 
we have 

(o_w o\, —2 =n... ‘2... ws} a2 

2 —W aa(t=— 2Qq, = = \at G+2Qq, = By’ 
(o_w\, — oP Ode, Oy , Oe _ 2.3 
ia ae" a. “= 2 ee oe 


From these equations it was shown by Squire (6) that (q,,¢,,q,) may be 
eliminated and a single equation for P obtained, viz. 


ed 





5- we a V#P +409, =0. (2.4) 

In discussing the motion in z > 0 let us assume that P and as many of 
its derivatives as necessary are given at z = 0 for all ¢ > 0 and that they 
tend to zero as x?+-y?-» oo. Further define P,, the Laplace transform of 
P with respect to t, by 


P(x, y, 2,8) = | e™ P(x, y,z, t) dt, (2.5) 
0 
1 c+ia 
; ae st 
so that P= a | P, e* ds, c> 90, (2.6) 
c—iso 


the path of integration lying to the right of all singularities of P,. If ulti- 
mately the flow is steady the value of P then is given by the residue of the 
pole of P, at s = 0. The equation satisfied by P, is, from (2.3) and (2.4), 


(+— w: =) VIP, +40) = = 40{%) , (2.7) 
t=0 

the right-hand side being the value of a(V2P) et when ¢ = 0: the initial 

value of V?P is zero because the relative motion is initially irrotational 

with potential ¢). On writing P, = P,—dpo, P, satisfies 


( ats =) VP, +4023 = 0. (2.8) 
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Since P, > 0 as x?+-y? + oo and is bounded in z > 0 we may write 
P, = ( eit dl ( ety Q)(l, m,z, 8) dm, (2.9) 
— 2 — x 


where /, m are real and Q satisfies 


(: = =) (ast o2Q) +40 F8 . 


. —=Q0 2.10 
dz} \ dz* dz? nan 


where o? = /?-++-m?, 
This equation is an ordinary linear equation of the fourth order and has 
four independent solutions, each of the form exp az, where a is a root of the 


quartic (s—aW)*(a2?—o?)+4Q2%2a2 = 0. (2.11) 
When s is large the four roots of (2.11) are 
+o,  (s/W)+ik, (2.12) 


k being equal to 20/W. 
Hence there is only one acceptable root a, in z > 0 leading to 

Q = F(s,l,m)e* (2.13) 
where F is an unknown function and a, is that root which tends to —o 
as soo. The other three roots can only be used to specify the solution 
in z < 0, and two of them correspond to waves carried along with the 
stream. To discuss the ultimate steady motion the behaviour of a, near 
s = 0 must be determined. Before doing this, however, we shall prove 
that if the real part of s is positive, the real part of a, is negative, and that 
the branch points of «9, gua function of s, are either on the imaginary axis 
or on the negative part of the real axis. 

First, «9 is a continuous function of s and has a negative real part when s 
has a large positive real part. Further, if the real part of a, is zero, we may 
set ay = 18, where f is real, 
ie = SF, 

Bo 
and so s is wholly imaginary. Hence a, can only cross the imaginary axis 
with s and, therefore, if the rea! part of s is positive, the real part of a, is 
negative. In addition, there cannot be any interchange between a, and the 
other roots in this domain of s. 

Secondly, the branch points of «, gua function of s, occur when two of the 
roots of (2.11) are equal. It may be shown after some elimination that this 


occurs when s = Wot). 2.14) 
At this value of s the roots are 

ol(oi—ki)!, ot(oi—kiyi 
and (of —k4)t{ —ki+ (of —otkt + kf)t}. (2.15) 
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We can make a, a single-valued function of s by requiring that a, -» —o 
as oo along the positive real axis and cutting the s-plane from the 
relevant branch points in (2.15) to infinity along straight lines parallel to 
the negative real axis of s. If o > k there are two branch points at real 
values of s, but from (2.15) the positive value of s corresponds to two of the 
roots with positive real parts being equal and is not a branch point of ap. 
Therefore, if « > k the only branch point of «9, gua function of s, is at 

s = —W(ot—k)}. (2.16) 
If o < k, however, there are two on the imaginary axis of s. 

From these two results we deduce that it is permissible to determine the 
behaviour of a) as s > 0 by considering positive real values of s only. Then 
the left-hand side of (2.11) is negative when a = 0 and positive when 
a = —0o so that there is always a negative real root which must be ay. 
In particular when s is small but positive the root is 

—(o?—k*}!+O(s) if o>k 

and is a gt Os?) if o<k. (2.17) 
Hence in the limit s > 0, a, > —(o?—k*)! if o > k and a > 0 if ao < &k. 
If in the solution for Q, the value of F in (2.13) with ?2+m? < k* is not 
identically zero, then the ultimate steady motion must include a cylindrical 
component. Therefore it does not tend to the undisturbed state as z > -++-00 
and the premiss on which the hypothesis is based does not hold. On the 
other hand, if there is no cylindrical component, 


F=0 if o<k, 


and in the special case of axial symmetry we have immediately the form 
of the hypothesis in (1.6). 

The same argument may be applied to regions far downstream of the body 
by considering the three roots of (2.11) with positive real parts when s has a 
positive real part. It may be shown that as s > 0 through positive values 
these roots tend to 0,0, (2—k)! if o>k 


and 0, ti(k?—o*}t if ao<k. (2.18) 


Thus, even if ultimately there is no two-dimensional component at large 
distances downstream there are still wave-like components coming from 
the conjugate imaginary roots when o < k. 

Finally, it is noted that the branch points are an indication of any 
inherent tendency of the solution to become unstable for large ¢t. A branch 
point with a positive real part will make a contribution to P which increases 
exponentially in amplitude while the contribution from a branch point 
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with a negative real part will tend to zero exponentially with time. Branch 
points which are wholly imaginary correspond to oscillatory motions 
which may increase or decrease algebraically in amplitude. The solution 
of the present problem when k = oo is already known (4): in it the branch 
points all occur on the imaginary axis in agreement with (2.16) and, further, 
the contributions to the solution from them.are algebraically small when t 
is large. It is likely therefore that at finite values of k the imaginary branch 
points make small contributions to the solution at large times. The con- 
tribution from the real branch point is exponentially small. There is a 
possibility that F may have a pole or branch points at a different point 
from those discussed above but again arguing from the solution when 
k = o, this is not likely. It is concluded therefore that the solution ahead 
of the body is ultimately steady. 


3. The stream functions which satisfy (1.6) 
We have from (1.5) 


,y=r ( A(c)J,(or)exp{—z(o?—k*) do if z>0 (3.1) 
k 


to satisfy (1.6). When r is small 


b = yr* [ cA(o)exp{—z(o*—k*)'} do = fr [ B(r)e-* dr, (3.2) 
k i} 
where 7? = o?—k? and rA(c) = B(r). 
Finally, expanding B in a power series about tr = 0 we have 
yp = kr? ¥ C,z-* 
in which the s are arbitrary positive numbers and C, arbitrary constants. 
Or, if (R, @) is a point in the fluid ahead of the body where 
Reosé = z, Rsin#é = r, (3.3) 
then yy = SC, R®-* when @ is sufficiently small. (3.4) 
It is required to find a solution of (1.2) such that (3.4) is satisfied and 
such that, as @ > 7, 4 > 0, so that the streamlines close up behind the body. 
Now the general solution of (1.2) which is bounded at both @ = Oand 6 = 7 
is a sum of terms of the form 
(1—p*) RV, (AR) P)(u) and (1—p?)RAJ_,, _(kR)P,(u), (3.5) 
where » = cos 6, n is a positive integer, J, ,, is a Bessel function, and P, (2) 
is a Legendre polynomial. The solution corresponding to n = 0 is 
RJ, (kR)(A+ Bcos 8), (3.6) 


where A and B are constants, and only vanishes at » = +1ifA = B= 0. 
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We now express that term of (3.4) which is proportional to R‘++, i < 3, 
together with any other terms of (3.4), of smaller order as R > oo, as may 
be convenient, in terms of Bessel functions. The series of Bessel functions 
obtained may then be generalized into a solution of (1.2), which satisfies 
(1.6), on making use of (3.5). Consider 


t . 
ba _ (4k R)s+t( — 1 s 
ig = Ri(1—p?) Ps eta (3.7) 





Provided only that 2t+-¢ < 3, ¢;, has the form of (3.4). Using the known 


‘denti 
we — (m+ 2u)(m+u—1)! 


u! 


(kR)™ = 





Iin+2u(kR) (3.8) 


u=0 
and mathematical induction, it may be proved that 
< (i+2u)(it+u+t)!(—1)Jj.0,(kR) 
= Ril—pv2JA(LR 1—y?2) Ri 2u 
(3.9) 
Hence the solution of (1.2), which behaves like ¢;, when @ is small, is 
RAL — p(k R) Piao) 
(¢+4)(¢—4) 





u=t+1 








+ Ri(1—p?) > (i-+-2u)(t+-u+t)\(—1)Jj 2k R) Pi 2u-4(#) 
# ui -u)(u—t— 1) (i+)! 8+ 2u4+-4)(+2u—4) 


(3.10) 


From the properties of P,() we know that n must be an integer and from 
(3.6) i+ 2u—4 cannot vanish for any allowable i, wu. Although 2¢ need only 
be less than §—1, it is convenient to take it as large as possible so that a 
Bessel function of negative order occurs only in the first term of the solution. 
First take therefore 


a 
u=t+1 


i= —}—2m, t=m and u=m+1+8 


where m, s-++1 are non-zero positive integers. Then (3.10) reduces, after 
multiplication by a numerical factor, to 


bom(R, 0) = Ri(1—p?)Pem(4)I4-2m(k-R)+ 


| 2 (48-+3)(2s)!(2m+1)! Ps,.1(u)Joq.4(kR) 
at i od > 8!(8+ 1)! m!(m—1)!(m+<s8-+ 1)(28+ 1 — 2m) 2+ 2m+1° 


(3.11) 





s=0 


Secondly, take i=}3—2m, t=m, u=m+1+s8 
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in (3.10). Then corresponding to (3.11) we have 
fom(R, 6) = RX — p*) Pom (HS y_om(kR)— 





~~ (48-+-5)(28-+-1)!(2m)! Po,(u)J) .26(4R) 
= 
aii 2, s!(s+1)! m!(m— 1)!(m+e8-+ 1)(2s— 2m +-3)228+2m+1° 
(3.12) 


The special case of (3.12) when m = | has also been obtained by Fraenkel 
(5) in a different way. 

A partial check on these formulae may be made by evaluating (3.11) 
and (3.12) when & is large. It is then found, for example, that the leading 
term of the series in (3.12) is 


2\'. 
(= ym RAL a p*)sgn BP oy, _4(#) +p} sn kR 
akR 
and cancels with the first term if 1 > 0. A similar result holds for (3.11). 


4. Flow past a sphere 

In this section we make use of the family of stream functions y,, satisfying 
(1.6) to determine the steady flow round a fixed sphere of radius a with its 
centre on the axis of rotation when at an infinite distance upstream the 
fluid has an axial component of velocity —W. Thus as z > « 


Y > —4Wr'. (4.1) 


From the previous section we know that the general solution of (1.2) 

satisfying (1.6) is 

Wai 
se 


= 


Y —4Wr? + b( R,O)A; (4.2) 
2 

and the A; have to be determined so that ¥ = 0 on R = a. Let us write 
(3.11) and (3.12) in the form 

py, = (1 —p*)Pi(u) RA, (KR) + PD» OK, 5( 1 —p*)Pi(u) Rt; (RR), (4.3) 

it 

where for convenience the numerical coefficients of the terms in the series 
have been replaced by a,;, half of which are zero. We must choose the A; 
so that 


(1—p*) > > {A,8;;J4_;(ha)+ A; a5} ,,(ka)}Pi(m) = (1—p?), (4.4) 
t=1j=1 


where 5,; = 1 if¢ = j and vanishes otherwise. Inverting the order of sum- 
mation, which may be justified a posteriori, 


> (1—p?)* Pi() Z {A;8;;J_4_(ha)+ A; a5; ,,(ka)} = (1—p*). (4.5) 


j=1 


Now the set of orthogonal functions (1—?)Pi(u), 7 = 1, 2, 3,... is complete 
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in the range |x| < 1 for all piecewise continuous functions which behave 
like (1—,*)! near z? = 1. Hence in (4.5) the A; must satisfy 


and the solution of this infinite set of linear equations determines ¥’. It has 


been solved numerically for ka = 1, 2, 3, 4, and the results are displayed 
in Table 1. 














TABLE 1 
ka 1, rer tae Se xh * | ties Bee 
I 0°903 —o'ol O13 
2°08 —o35 «| —oo! | 14 
3°81 —248 | -034 =| —O15 14 
4 —382 | —1092 | —4-36 ; = i.) oe to 
When ka = 5, A, = —15, but between ka = 4 and ka = 5 it changes 
sign twice since it is positive when ka = 4-49 where J,(ka) = 0. Thereafter 
it rapidly decreases, being —23 at ka = 5-2, —83 at ka = 5-7, infinitely close 


to ka = 5-76, and +31 at ka = 6. At ka = 5-76 the determinant of the 
coefficients of (4.6) when divided by the diagonal terms, which dominate it 
at sufficiently large values of i,j, vanishes without any of these terms being 
infinite. Although from a mathematical standpoint solutions can be 
obtained at higher values of ka it is physically inconceivable that this 
solution could oceur when ka ~ 5-76 and therefore also at higher values the 
mathematical solution must be rejected. Tt is concluded therefore that 
the problem is not properly posed if ka > 5-76 and, indeed, it would be 
surprising if it did not cease to be properly posed at a lower value of ka. In 
the discussion of the drag given below we shall see that the drag on the 
sphere is out of all proportion even at values of ka as low as 3. It is interesting 
to note that according to Taylor’s experiments (1) the character of the flow 
changed to include a cylindrical component when ka ~ 27. Although this re- 
sult was only obtained approximately so that it is not significantly different 
from 5-76 the agreement must be regarded as fortuitous since the physical 
breakdown may well occur before then. 
From Squire’s work (6) the pressure p in the fluid is given by 
k*y? 

p = const.—4}p(u?+ w*)—}p bad 4Q%pr?, (4.7) 
where p is the density. Hence since ¢ = 4Wa’sin?# and ‘Y = 0 on the 
sphere, the only part p, of the pressure which contributes to the drag D is 


p (er)? 
si £1) (4.8) 


5092.41 E 
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which comes from the longitudinal component of velocity. Hence 


7 +1 
ane 
D= [ Pp cos 8 2na*sin 6 dé = — | HS) du. (4.9) 
b 1—p*\oR 
0 -1 


Now from (4.2), using (4.3) and (4.6), 


l oy _ 3 ka Ji(ka)\ py (—)*14, 
yer” hgaliaeiealt (es Jest Pi(u +> Ti-(eay * 





= —p (—)'*1B,; Pi(u) say. (4.10) 


Hence, substituting in (4.9) and writing 
D = pna*W*C;(ka), 


where C, is the drag coefficient, we have 


cs +1 
Cp = — pp) J B, By —) p(l—p?)Pi(u) Pi(u) dp 
0 2 +1 
=D BB [ny P, Pi 
a BA 
i(i+1)(t+2) 
ee 4.11 
5 ate 


If ka is so small that terms O(ka)!* may be neglected, C,, may be calculated 
from A, and A, only. It is then found that 
A, = [J_,(ka)|",  B, = wkaA,J,(ka) and B, = 3A,[J_,(ka)]|" 
(4.12) 
(ka)* 7 : ' 
[1—(ka)?+-(ka)*—§(ka)®]+ O(ka)®. (4.13) 


so that Cp = ri 


For larger values of ka, the values of C, are given in Table 1 and may 
be compared with the value 

Cp = 8ka/3zn, (4.14) 

valid when ka is large, obtained in (4). In that solution (1.6) is inapplicable 

and the flow is wholly cylindrical. From a comparison between Table 1 

and (4.14) it is inferred that the values of C, obtained when ka = 3 and 


ka = 4 are too large and that in practice there would probably be a cylindri- 
cal component to the flow. 


The author is grateful to Mr. L. E. Fraenkel for his constructive criticism 
of the first draft of this paper. 
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HEAT TRANSFER THROUGH THE LAMINAR 
BOUNDARY LAYER ON A CIRCULAR CYLINDER 
IN AXIAL INCOMPRESSIBLE FLOW 
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SUMMARY 


This paper presents a method of calculating the distribution of rate of heat transfer 
into a laminar incompressible boundary layer from the exterior surface of a long thin 
circular cylinder, when the surface of the cylinder is maintained at a constant tem- 
perature and the flow is parallel to the cylinder axis; the temperature difference 
between the surface and the main stream is taken to be small enough to neglect 
buoyancy effects. A series solution, valid for small downstream distances from the 
nose, has been obtained already by Seban, Bond, and Kelly. This is now extended 
by deriving an asymptotic series solution, valid at large downstream distances, and 
bridging the gap between these two series solutions by an approximate solution, 
based on the method used recently by Davies and Bourne to calculate heat transfer 
from a flat plate. The calculation is used to demonstrate the effect of curvature 
and of Prandtl number on the local rate of heat transfer at various downstream 
distances by comparing with the corresponding flat plate results. 


1. Introduction 


In order to estimate the influence of curvature on the calculated values 
of local rate of heat transfer from heated surfaces into a given flow, the 
problem of forced convection from a heated circular cylinder into the 
surrounding axisymmetric, incompressible, laminar, boundary layer pro- 
duced by a uniform stream was considered first by Seban and Bond (1). 
The front of the cylinder is supposed to be smooth enough to prevent 
boundary layer separation, and the variation with downstream distance 
of the local rate of heat transfer Q from the external surface of the cylinder 
is investigated. Seban and Bond give an exact solution of the appropriate 
temperature equation in the form of a series, valid for small values of 


ve/Ua? where U is the uniform mainstream velocity, x is downstream 


distance, v kinematic viscosity, and a the cylinder radius. Important 
numerical corrections of this work by Seban and Bond were later given 
by Kelly (2), and we shall refer to their solution as the Seban—Bond—Kelly 
solution. The leading term of this series solution is equivalent to the flat 
plate solution, while the succeeding terms represent corrections for the 
effect of curvature which becomes increasingly important at large down- 
stream distances. The Seban—Bond—Kelly solution extends the flat plate 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958] 
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solution, which is valid for vglues of va/Ua? of less than about 0-0001, up 
to approximately va/Ua? = 0-04. 

We first obtain an alternative solution for the region vz/Ua? < 0-0001 
by utilizing the Pohlhausen type of velocity profile given by Glauert and 
Lighthill (3) and extending a method of treating the heated flat plate 
problem described by Davies and Bourne (4). In the case of constant 
surface temperature the numerical result is within 2 per cent of the exact 
flat plate value and suggests that the Glauert—Lighthill profile is a suitable 
approximation to the basic flow distribution in this context. 

We next obtain an asymptotic series solution, which is valid at large 
values of vx/U’a®, based on the analysis described by Glauert and Lighthill 
(3). The series is expressed in terms of inverse powers of 8 = In(4va/Ua*) 
up to the term involving 1/8%. As in the skin friction problem (3) this is 
sufficient to calculate Q for va/Ua? > 100, although the results given by 
Glauert and Lighthill suggest that the series can be used, provided the 
Prandtl number is approximately equal to one, up to about vx/Ua* = 10. 

In order to bridge the gap between the Seban—Bond—Kelly solution and 
the asymptotic solution an approximate method is required. A Karman— 
Pohlhausen approximation could be employed, following the method given 
by Glauert and Lighthill, and the distribution of C computed by numerical 
integration of a first-order linear differential equation for each prescribed 
Prandtl number. However, this method of approach is only applicable in 
the case of constant cylinder temperature. We use instead an alternative 
method which is applicable to an arbitrary distribution of cylinder tem- 
perature, and, in the case of constant cylinder temperature, it leads to a 
fairly simple computation, involving the numerical evaluation of quadra- 
tures at appropriate positions on the bridge. We begin with power law 
approximations of similar type to those used by Davies and Bourne, 
involving deviations from the Pohlhausen profile of similar magnitude to 
those experienced in the flat plate problem (4). The ensuing temperature 
equation is transformed into a von Mises form and the appropriate boundary 
conditions specified by regarding the heated cylinder as an assembly of 
continuous uniform circular sources, each source being positioned normally 
to the cylinder axis and of radius a. A solution of the temperature equation 
is then obtained for one circular source and the distribution of source 
strength (or local rate of heat transfer) is determined by solving an appro- 
priate integral equation, so that the prescribed cylinder temperature is 
achieved. We find that the calculated value of Q at va/Ua? = 0-04 is 
within about 3 per cent of the value given by the exact Seban—Bond-Kelly 
solution, and the values of Q for 0-04 < va/Ua* < 10 continue into those 
computed from the asymptotic solution. For reasons discussed in section 5 
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it is probable that the values of Q based on this approximate method are 
in error by no more than about 5 per cent. 

The influence of curvature is displayed by comparing the calculated 
numerical values of Q, for air, with those given by the flat plate solution, 
and the dependence of Q on o, the Prandtl number, is also discussed. The 
thickness of the boundary layer, at small values of va/Ua?, is very small 
compared with the cylinder radius, and consequently Q varies with o* as 
in the flat plate case, but as va/Ua? increases the value of the index of o 
decreases from }, and at very large values of va/Ua?, Q is independent of c. 
When the numerical value of o is about 1-0, the variation with va/Ua? of 
the quantities given by (a) [local skin-friction on the cylinder/local skin- 
friction on the flat plate], and (5) [local rate of heat transfer on the cylinder/ 
local rate of heat transfer on the flat plate] are similar, but when o differs 
markedly from 1-0 this is not true. Thus at very large x-values, the asymp- 
totic solution shows that the ratio of (a) to (b), at a prescribed z, is given 
by o-*, while at small values of vx/Ua? it is 1-0. 


2. Solutions for small downstream distances 
The boundary layer equations are 
ou cu vof{ ou 
6 — + 8 — = - —1 fF — (1) 
Ox er ror\ or 


for the velocity distribution, 


oT oT «ela 
Oo 0 ae — Pe 
Cx or ror\ or 


for the temperature distribution (neglecting frictional heating), and the 
continuity equation is 


© (rv) 4 © (ru) = 0, (3) 
cr Ox 


where x denotes downstream distance, r the radial distance from the axis, 
u and v the downstream and radial velocity components respectively, 
7 the temperature, and «x the thermal diffusivity. 

The Seban—Bond—Kelly series solution for heat transfer, referred to in 
the Introduction, is valid for va/Ua? < 0-04 and is given, in the case of 
constant cylinder temperature, in the form 

Q Ua*\}/ va \} | 
——* ___ = 0-664mo#{ ——}®/ 1 +-2-30{-— }* — ..., 4 
WET) = Uai) ~""| (*) 
where Q is the local rate of heat transfer (per unit length of the cylinder), 
o is the Prandtl number, k the thermal conductivity, 7, the cylinder 








Pa Tansee atin eesti 
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temperature, and 7, the air stream temperature. The leading term in 
(4) represents the flat plate solution for a Blasius layer. 

We now consider an alternative method of approach which suggests the 
approximate treatment used for larger values of vz/Ua*. We adopt the 
Pohlhausen profile used successfully by Glauert and Lighthill (3). The 
order of magnitude of the error involved, due to the use of this profile, 
in the heat problem can also be estimated from this by comparing the 
ensuing solution with the exact solution (4). The Pohlhausen profile is 

u 1 
U (2) 


and i: =1, fory >3, (6) 


where 6 is the ‘boundary layer thickness’ which appears in Pohlhausen 
treatments and y denotes distance from the cylinder surface. The parameter 
. is given numerically as a function of x by Glauert and Lighthill. We 
detine the stream function %, to satisfy (3), by the relations 

a =ru and op = —rv. (7) 
ér ex 


In(l+y/a), for y <8 = a(e*—1), (5) 


Expanding the logarithmic form in (5) and retaining only the first term, 
the velocity profile, at very small va/Ua* values, is given by 


~ o where { = %, (8) 
U a a 
and, using (7), the stream function ~% becomes 
Ua? 
op my cae - ¢?, (9) 
aX 


The von Mises method (5, p. 126), used in the case of the flat plate (4), 
is easily extended to the cylinder aga and the basic equation (2) simpli- 


ties to the form P 
ox a} 
Neglecting squares and higher powers fe { we have, using (8), 
ru = a*Ua-l, (11) 
and hence, using (9) and (11), equation (10) now becomes 
eT _ = (24U tata) — SS. (12) 
oa 


Using the Seban-Bond—Kelly solution, it can be shown that, when 
ve Uat*+>0O, 


o Ua? 


Bs o-ss2( 72) 4, (13) 
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and on substituting X = z?, (12) reduces to 


eT a/,,2 
tS yi et 14 
ax ex exp | sine 


where c = 243-1(0-332)'Utalav. 

Equation (14) is identical in form with basic equation (11) discussed by 
Davies and Bourne (4), and their analysis applies with Q/27a replacing Q, 
the rate of heat transfer per unit length of a line source. We obtain 

RET = 0-678r0}(" =) (15) 
This result is within about 2 per cent of the exact flat plate solution. The 
small error is due to the logarithmic approximation of the Pohlhausen 
profile but, as discussed by Glauert and Lighthill, the error on this account 
is likely to decrease as vx/Ua? increases; for this reason the logarithmic 
profile (5) is used in the following sections, for all values of vz/Ua?. Solution 
(15) is only valid down to va/Ua* = 0-0001, when higher terms in ¢ 
become significant, so that an alternative method must be used. 


3. The asymptotic series solution for large downstream distances 

Following the method of treating the velocity boundary layer equation 
(1) given by Glauert and Lighthill (3),+ we now obtain a solution of the 
temperature equation (2), which will be valid at large downstream dis- 
tances where the boundary layer thickness is much larger than the cylinder 
radius a. Glauert and Lighthill consider a series expansion for the stream 
function 


AME) AE), | , 

= vaf(é,2) ~ vel fol) + ff) +. (16) 

where & == a and B= inf) (17) 
so that u=i3Uf’ and v= -(f'—f—xf.); (18) 


primes and x suffixes denote partial differentiation with respect to € and x 
respectively. They show that 


fo = 2, (19) 
f, = 2€ Ei(—£)+2e*—2, (20) 


f2 = 2e-* Ei(—é)—4 Ei(— 2¢)—{Ei(—£)}?+-4 Ei(—é) In €+ 
+6 Ei(—£)—6El(—€), (21) 


é é 
where  Ei(—€) = [ e“t4dt and El(—£) = | Ei(—t)t— dt. 


i 2) 


¢ An asymptotic series solution for skin-friction has also been obtained by Stewartson (6). 

















HEAT TRANSFER THROUGH THE LAMINAR BOUNDARY LAYER 57 


The similarity in form of the velocity and temperature equations 
immediately suggests that an asymptotic series solution, of the same nature 
as the Glauert—Lighthill solution for the velocity profile, exists for the 
temperature profile. Accordingly we consider the series expansion 


Tat, = 62) = ge) +29) 4 4... (22) 

the temperature boundary conditions being expressed in the form 
g>l asé>o, (23) 
and g>0 asé->e%, (24) 


Using (18) and (22) the temperature equation (2) becomes 
af '9,—f9' —xf,g' = 20-19 + 20-%f9”. (25) 


Substituting the series (16) and (22) and equating inverse powers of 8 we 
obtain the equations 


E9o+Go+ bofogo = 9, (26) 
Egitgit4ogohitgifo) = 9, (27) 
£93+92+40(fego—fidot+tomt+hgithog2) = 9, (28) 


and in general 
n n—1 
EGn T In T bo| > Sm Gn- m— 5, Sn Fn-m—1—ImJ n-m-1)} = 0. (29) 
Following closely the analysis of Glauert and Lighthill for the velocity 


profile, we now transfer the inner boundary condition (24) to € = 0, since 
es Gas-as are independent of 8. Now it is easily seen that near ¢ = 0 
equation (25) has the form 

gg" +9 = 0, 
with solution g = C+ Diné. 


Hence, near £ = 0, each separate term in the series (22) must have the 
form 


g~c,+d, In, (30) 


where c,, and d,, are constants. Using these forms in (22) with the condition 


(24) we obtain a 
C,—Pa, 
>, ar =o (31) 


n=0 
so that d, = 0, d, = Ce, dy = ),..., d,, = Cy_+. (32) 
The inner boundary condition is thus expressed by (30) and (32) as € + 0. 
The condition at infinity (23) becomes 
Jo > 1, Jn>0(n>0), as&>o. (33) 
Substituting (19) into (26), we find the equation for g, is 
£9 +Got+obGo = 9, 
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the solution of which, under the conditions g, > c, as > 0 and g, > | as 
£ -» 0, is seen by inspection to be 
Jo = 1. (34) 
This represents an isothermal solution 7’ = 7, and corresponds, as we 
would expect, to the solution f, = 2 which represents a uniform stream 
of velocity U in the velocity case. 
Equation (27) now becomes 
gi +g, +08, = 0, (35) 
with associated boundary conditions 
g,~C,+Iné as€>0, and g,>0 asé>o—m. (36) 
The appropriate solution of (35) is 
g, = GEi(—e6), 
where G@ is an arbitrary constant. As £— 0, it may be shown that 
Ei(—of) ~ Iné+(y+lIno). 


Hence, using (36), we find that 


c¢, = y+lIno, (37) 
and g, = Ei(—e6). (38) 
Using (19)-(21), (34), and (38), equation (28) becomes 
Eg -+9+otgs+{o Bi(—E)e-* + 0€-1e-1 + +o Ei(—o£) —of-4e-*F} = 0, 


with the boundary conditions (39) 


Jo™ Co+(y+Ino)Iné asé>0, and g,>0 as&>o. (40) 
Integration of (39) yields 
Jo = e-% Ei(—£é)—(1+0) Eif—(1+o)€}+(1+ 0) Ei(—o€) In — 
c 
—(2+-0) El(—o€)+ H Ei(—s€)—o | e-*Ei(—t)t-"" dt, (41) 
where H is an arbitrary constant. It is easily shown that each member 
of (41) tends to zero as € tends to infinity, and the boundary condition 
as € + 00 is satisfied. In order to apply the inner boundary condition we 
use the result, quoted by Glauert and Lighthill (3), that as é > 0 
El(—£€) ~ 4(y+1n €)? +42? (42) 
and the result, obtained in the Appendix of this paper that, as € > 0, 


é 
| e-% Ei(—t)t-! dt ~ }(y+1n €)?+47°—S,, (43) 


: 3 
here GEE ae ae 44 
whe 2 ata (44) 


Co 
» 
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The series S, is convergent and therefore directly summable for |o| < 1 
which is the case considered in section 5 of this paper (for air o = 0-72). 
If o > 1, the series diverges but is nevertheless summable by Euler’s 
method (see, for example, Hardy (7)).t It follows that, as ¢ > 0, 


jg, ~ {H(y+lno)—yo—y*({1+oe)—}(1+oe)n?—(1+¢) In(1+o)— 
~}(2+0)(2y+Ino)Ino+oS,}+{H—o—(1+e)y—Ino}ln€, (45) 
and applying the inner boundary condition (40) we fiad 


H = o(1+y)+2(y+lno), (46) 


d, = y?+(e0+2y)Ine+(1—4o)(Ino)?— 
—(1+e)In(1+0)—}(1+0)z*+-o8,. (47) 
When o = |, these reduce to the equivalent results obtained by Glauert 
and Lighthill (3). 

Proceeding in this way accuracy to any inverse power of 8 can be achieved 
but the integrations clearly become extremely complicated. In the case 
of air we find, following Glauert and Lighthill, that the results obtained 
up to this stage are sufficient to compute the rate of heat transfer for 
ve L'a? > 100. In the skin friction case the frictional drag +r, calculated 
with three terms in the series for f, was still in agreement at va/Ua* = 10 
with the result suggested by the Pohlhausen method. It therefore appears 
probable that these first three terms are sufficient to estimate Q for 
ve Ua® > 10 for the flow of air past a heated cylinder. 


The local rate of heat transfer, per unit length of the cylinder, is now 
given by 


¢ ae mere 
gen = ef) = tears ~ eel 


n=0 


sé 0, ég’, > d,,, and hence 


Q Sd 
a eae _ 
k(T,—T) "Le 


The calculated values of d,, d,, d, show finally that 


_(y+ino) , 


wn Gat «<b. Siewtaaee 


at B* 


a+ 2y)lIno+(1—4$oe)(In o)?—(1+-¢) In(1+oe)—j(1 tome) 
p 





_ 


+0(34) (50) 


+ S, is tabulated in the range —1 < o < 5 (see E. O. Powell, Phil. Mag. 34 (1943), 600). 
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It is also of considerable interest to note that at sufficiently large values 


of x (or 8), Q is independent of oc, so that, using (50) and (4), the value of 
the ratio 


F _ Q at downstream distance x for a cylinder 4nB-! t 
= ene aeennEne ~ — - ——_—_——____—_———, 
1 Q at downstream distance z for a flat plate — 0-66470!(Ua®/vz)! 


4 


the value of the corresponding quantity of skin-friction, using the 
Glauert—Lighthill asymptotic solution and the Blasius flat plate solution, is 


__ 7 at downstream distance 2 for a cylinder 4nB-! 
> 7 at downstream distance «x for a flat plate 0-6642(Ua?/vax)* 


Hence, for large values of vx/Ua?, 


while at very small values of vz/Ua?, this is, of course, unity. These results 
are of importance in calculating particular cases, as discussed in section 5. 


4. Anapproximate solution for intermediate downstream distances 


Lighthill (8), in treating the heated flat plate problem, found that a 
linear approximation to the exact velocity profile in the neighbourhood of 
the surface led to a reasonably accurate result for local rate of heat transfer 
Q in certain cases of the velocity form U = cx”, but in other cases, when 
the exact velocity curves deviate considerably from the tangent at the 
origin, a considerable degree of error ensued in the calculated values of Q. 
However, this error was reduced to very small proportions by Davies and 
Bourne (4), using power law representations of the exact profiles. This 
suggests that an approximate method based on similar premisses is likely 
to yield sufficiently accurate results here and we show in section 5 that 
this is so. 

Using the continuity equation (7) and the Pohlhausen velocity profile 
(5) we first write ap aU 

é “ti x 


on 


where 7 = r/a. ee yields 


nin», (52) 


r 


$=5 — {(2In »—1)n?+ 1}, 


so that 4 = 0 when 7 = 1. The temperature equation (10) in von Mises 
form then becomes 
eT = xa®U ? f(y or" 
ex Par ike 8 


In order to proceed further and apply the method of sources, used 


(54) 
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successfully in the heated flat plate problem by Davies and Bourne (4), we 
“= (2In 7—1)y2+1 = X(y), —-o*Ing = F(n), (55) 


and suppose that a good power law representation of the variation of Y 
with X can be constructed in the form 


Y = BX, (56) 


where B and ¢ are constants. In calculating Q at a given downstream 
distance x, the numerical values of B and t are chosen to give a good fit 
over the range of » corresponding to the thickness of the boundary layer 
at this 2, value. This procedure is then repeated for other values of 2,; 
the numerical details and implications in the case of air flow over the 
cylinder are discussed in section 5. Substitution of (55) and (56) into 
(54) now gives 
oT . a 2 
— (4! Ba2@U Aq) gt : 57 
oP = wnt Sy) (51) 
Equation (57) next reduces to the form of the basic equation (11) of 
(4), if we write 
[ {a(x)}¢- de = X’, (58) 
0 


where a(x) is given numerically by the Glauert—Lighthill results. We now 


have : 
eT 2 /(,,@ 
pte ey eM Piel 59 
ex’ (va } (59) 


where c = 4! Ba%1-UG.«, 


Suppose now that acontinuous uniform circular source, of radius a, normal 
to the generators of the cylinder, is stationed on the surface of the cylinder, 
yy = 0, at x = x», and emits Q units of heat per unit of time. The analysis 
of the flat plate approximate solution (4) follows with Q replaced by Q/27a 
and the distribution of local rate of heat transfer may be calculated for 
any prescribed downstream distribution of surface temperature. The result, 
when the surface temperature of the cylinder is independent of z, is 


@ A224) gin| 7 1 \5 | 
- - = Ze ) eciteles 2—t (2-4 R1(2-4*) 
KT) sin( 7 )P(g =) 2-1 Bx 


2—t 


va|Uat ; 
a (ra a { e-0 ae aaa, (60) 
0 


Ua? 


\° 


where a’ = va/Ua?. 
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5. Calculated distribution of local rate of heat transfer in air 
To illustrate numerically the formal results obtained in the previous 
sections, we now consider the case of air flow over the cylinder, for which 
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Fic. 1. Distribution of Q, local rate of heat transfer, along the surface of a circular 
cylinder, heated to a constant temperature and positioned with its axis parallel to 
a uniform stream of air. Q, rate of heat transfer per unit length of cylinder; k, 
T,), temperature difference between cylinder and main 
main stream velocity ; 


thermal conductivity ; (7; 
stream ; v, kinematic viscosity ; 7, downstream distance ; U, 
a, radius of cylinder ; — — — — values of Q based on flat plate solution ; ———-— 
values of Q based on an exact representation of conditions of flow near the surface ; 
———— values of Q based on the Seban-Bond—Kelly solution (on the left) and the 


asymptotic solution (on the right); .......... values of Q based on the power 


law approximate solution 


we take o = 0-72. All the results discussed in this section are exhibited 
graphically in Fig. 1, together with the corresponding results for a flat plate. 

As in the skin-friction problem, the probable limit of reliability of the 
Seban—Bond—Kelly temperature solution is at a value of vx/Ua? of about 
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0-04. The limit of reliability of the asymptotic series solution is probably 
at about vz/Ua® = 100, where the third term in the series is about 9 per cent 
of the first two terms together (taking o = 0-72). Atva/Ua? = 10, itis about 
28 per cent of the first two taken together; these results are similar to those 
discussed by Glauert and Lighthill. However, since a is of the order of one 
for air and, since the Glauert—Lighthill asymptotic solution appears to be 
a good approximation to the true skin-friction variation up to about 
vx/Ua® = 10, it is probable that the asymptotic series solution (50) for the 
temperature problem is a good approximation to the true values of Q up 
to about va/Ua* = 10. 

The approximate solution described in section 4 is now used to bridge 
the gap, 0:04 < vz/Ua* < 10, between the results computed from the 
Seban-Bond—Kelly series and the asymptotic series. The computation of 
Q(x) by this approximate method involves (i) the evaluation of the para- 


vr/Ua* 
meters B and ¢, and (ii) the numerical evaluation of | a - dz’, using 


0 
the distribution of a given numerically by Glauert and Lighthill (3), and 
Q has been evaluated at several points on this ‘bridge’, the same procedure 
of power law fitting being used in each case. At a given intermediate value 
of vx/Ua?, the appropriate a value is noted from the Glauert—Lighthill 
results, and the thickness of the momentym boundary layer 5 evaluated 
using (5). The thickness of the temperature boundary layer 4,,, fora = 0-72, 
is then given approximately by the boundary layer equations (1) and (2) to 
be 6/(0-72)t. The values of X(n) and Y(»), introduced in section 4, are 
computed at intervals up to n, = 1+8/a, and in the uniform stream region 
above », the values of X(n;)+2a(n?— 1%) and 7*a, which correspond to 
X and Y respectively when y > 3, are computed from » = 1+-8/a to 
7 = 1+6,,/a. Values of B and tin the power law Y = BX‘ are now chosen 
to obtain zero deviation from the exact Y = f (X) relation at the mid-point 
of 5,, and to minimize the deviations elsewhere. In this way the best fit 
is obtained over the central bulk of the boundary layer (where most of the 
flux of heat takes place). The deviations are small, being 4 per cent at 
most over 95 per cent of the range at the smaller intermediate values and 
10 per cent at the larger values. Numerical values of B and t appropriate 
to a given vx/Ua? in the intermediate range do not of course yield a good 
approximate description of Y = f(X) at points far upstream. However, 
the calculated value of Q at a chosen distance z is likely to depend critically 
on the distribution of flow properties in the neighbourhood of x and is not 
so dependent on conditions far upstream, the main contribution to the 
temperature near the surface coming from the sources in the immediate 
upstream neighbourhood of x. The discontinuity in the slope, éu/éy, 
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associated with the profile (5) at y = 5 may also be considered an objection- 
able feature, but, as noted by Glauert and Lighthill, this does not necessarily 
lead to serious error in Pohlhausen treatments. The error is likely to be 
small for values of o only slightly less than 1-0 (as in the case of air, con- 
sidered in this section), since most of the heat flux across a plane normal 
to the surface will take place in the region y < 5. In the case of fluids, 
whose Prandtl numbers are in excess of unity, this difficulty does not enter 
at all, as 6,, < 5. 

We note that this approximate treatment, which has the advantage of 
involving fairly simple computation, is shown to lead to reasonable numeri- 
cal results (probably to within 5 per cent of the true values) for the following 
reasons: 

(i) the computed values of Q, based on this power law method, are found 
to be within 3 per cent of the Seban—Bond-Kelly results for the smaller 
values of va/Ua?; 

(ii) these values of Q run smoothly into the numerical results given by 
the asymptotic series solution at very large values of va/Ua?; 

(iii) all the values computed in the intermediate range of vx/Ua?* are 
found to be within 5 per cent of the results suggested by considering the 
ratio R,,/ Rp defined in section 3. The analysis described in section 3 shows 
that the numerical results for the ratio R,,/ Rp are likely to vary between 
1-0, at small va/U’a® values, up to 1-12 (for o = 0-72), at very large vx/Ua? 
values. At va/Ua? = 100 the asymptotic series solution for heat and skin 
friction shows that R,,/ Rp, = 1-06, so that over the intermediate range of 
va/Ua® the computed values of R,,/R, should lie between about 1-01 and 
1-05. This serves as a useful check on the numerical results given by our 
approximate method for R,,/ R, and we find that they are all within about 
5 per cent of this small range of probable values. Consideration of (iii) above 
could, of course, be used to deduce approximate values of the interpolation 
function needed. The power law analysis, however, is capable of applica- 
tion to any prescribed downstream distribution of surface temperature, and 
by its use formulae describing the temperature distribution in the airflow 
can easily be derived; the dependence of Q on a is also conveniently deduced 
by its use. 

In the case of the flat plate problem Q depends on o! at all values of 
va/Ua*®. The effect of curvature, however, on this dependence is found to 
be considerable. At small vx/Ua? values Q on the cylinder is, of course, 
proportional to o as in the flat plate problem, but the index of o decreases 
as downstream distance increases. This index is given by (1—t)/(2—t) in 
the analysis of section 4 (equation (60)); e.g. at vz/Ua? = 0-003, 0-06, 0-3, 
and 1-6, (1—t)/(2—?) is given by the power law results, for air, to be 0-31, 
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0-27, 0-21, and 0-18, while at vz/Ua? = o it is zero. These numerical 
results are likely to be approximately true for fluids whose Prandtl numbers 
are near 0-72. 

Finally, it is of some interest to calculate Q(x) over the whole range of 
vx Ua® by the method of section 4, using values of the parameters B and t 
which ensure an exact representation of conditions near the surface, as in 
Lighthill's linear approximation for the flat plate problem. The increasingly 
logarithmic behaviour of the velocity profile in the cylinder case, as x 
increases, results in a rapid divergence from the tangent linear approxima- 
tion to the profile near the surface and very good results cannot be expected 
in this context. As r >a (» > 1), it is easily seen, using (55) and (56), that 
t—4and B—1/v2. When va/Ua? = 0-0001, the value of Q calculated in 
this way is 2 per cent greater than that given by the exact flat plate solution 
and is in agreement with our approximate solution. At vz/Ua? = 0-001, 
Q is in excess by 4 per cent of the exact flat plate value, whereas the Seban— 
Bond-Kelly solution (and also our approximate solution at this stage) leads 
to an excess of 7 per cent; this is also the result obtained in the skin friction 
case. The results based on this approximation are seen from Fig. 1 to be 
seriously in error at larger downstream distances. At va/Ua? = 100, the 
result, though an improvement on the flat plate solution, is nevertheless 
still considerably below the true value calculated from the asymptotic 
solution. 

The effect of curvature on Q is clearly seen in Fig. 1 and shows, for 
instance, that the numerical value of Q calculated from the asymptotic 
solution, at vx’ Ua? = 100, is ten times the flat plate value. 
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APPENDIX 
In section 3 we made use of the result that, as > 0, 
g 
| e~°t Ei(—t)t-! dt ~ Hy+In€)*+y4y7*— §, (61) 
oc oe o& 


where S, ot ge To derive (61) we first integrate the left-hand side by 

parts, F 

{ e-°t Bi(—eye1 dt = HEi(—£)}%e-° — (1a) | (Ei(—2)}%e- dt. (62) 
@o 


x 


5092.41 F 
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€ 
Hence as £ > 0, e~* Ei( —t)t-! dt ~ H(y+In€)*+ J, (63) 


> 6) 


@ 
where J = }(1- a) f {Ei(—t)}%e"--" de and may be shown to be convergent for 
0 


a > —1. Expanding the exponential term in J and integrating by parts yields 
© wo t 
= : - r+lgrt+le—t Bj —t) dt 
I = }(1—¢) { ar (Ei —t)}* de - -{ Dt . 
or oe 
Hence 
I= — f fett-ot—1e-1¢ *Ei(—t) dt 
0 
x 
- lim {e~°tt-1 Ei( — t)—e—*t-! Ei( —t)} de 
Xo » 
ro * 





r+ lyr A r+ler By 
i Mie { >{t {(—o anti t) (—1) +t" Ei(—2)| és. 
Xo +1)! (r+1)! J 


z0 z r=0 


y [ (—o)"*12" Ei( —t) dt (—o) (+) 
Now ee 
(r+1)! (r+1)? 








by repeated integration by parts. Hence we find that 


{ : — (—o)"t? rr |. 
r= >(- m 4-- a = 7a (64) 


and the required result (61) fitinte We note that the result (42), quoted by Glauert 
and Lighthill, follows from (61) as the special case o = 0. 
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HYPO-ELASTIC POTENTIALS 
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SUMMARY 
For some hypo-elastic materials it is possible to introduce a scalar potential 
analogous to the elastic potential or strain energy function of elasticity theory. 
Conditions necessary and sufficient for the existence of such a potential are derived. 


1. Hypo-e.asticrry, as defined by Truesdell (1, 2, 3), is the theory of 
materials with constitutive equations given by relations of the form 
bis = bgt Oia — be Mee tty tee = Asjem Lem 

in cartesian tensor notation. Here ¢,; is the stress tensor, v, the velocity 
vector, 24h, = UpmtUmer Aijem 8 @ tensor invariant of t;; such that 
A ijcm = Ajizxm = Aijmes 2nd the ‘dot’ denotes the material derivative. 
For our purposes, it is slightly more convenient to use an equivalent 
formulation which has been emphasized by Noll (4, section 14) and Thomas 
(5, 6), namely 

Dts 2 = by —tywy—ty on = Bijymd (1) 

ry ij ik ‘jk “ik tjkm “km? 
where 2w;; = v,;—v;, and B;,,,, is a tensor invariant of t;; with the same 
symmetry as A;;,,. We assume that B;,,,, is a continuously differentiable 
function of t,,, in some region R of stress values. Also, in R, B;,,,., is required 
to possess an inverse C;,,,,, which is a tensor invariant of stress with the 
same symmetry as B;,,,, such that 


2C; jkm Bimre = oes 2 Bistem© ‘kmrs ~~ = 8,55. +8;,5 ir’ (2) 
It exists if and only if Dt,;/dt = 0 implies d,,; = 0. It is to be understood 
that values of ¢,; not in R are henceforth excluded. 
Our purpose je to characterize the class of hypo-elastic materials for 
which there exist scalar invariants ¢ and y of t,; such that the equations 
tidy; = op (3) 
hold whenever ¢,; and v; ; satisfy (1). When ¢ and ¢ exist, we call y a hypo- 
elastic potential. Thomas (7, 8) has considered equations similar to but not 
identical with (3). 
Noll (4, section 15) has shown that hypo-elasticity includes as a special 
case the more familiar theory of elasticity for isotropic materials which 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958] 
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assumes that stress is a function of strain, derivable from a scalar strain 
energy function. For the materials described by both theories, (3) holds; 
¢ can be taken to be the mass density and y the strain energy per unit mass, 
these being expressed as functions of stress. 


2. Any scalar invariant of ¢;; is expressible as a function of the three 


invariants - 
I, = ty, I, = tastes, Ty = bij tin tes 
It is easily verified that 
Iy Dt;; 
ty = 24, =o (N = 1,2,3), 
at;; ét;; dt 
from which it follows that (3) can be written in the form 
Op y ~ Dt; ep 
t... dp = = — B.... Gen: 4 
km “km : hip = ¢ — ét;; ti; *9= —— =¢.- at; ijkm “km ( ) 


In order that+ (4), holds for all t;; and y consistent with (1), it is necessary 
and sufficient that ‘ 
=¢ ob RB. (5) 
nd 
To see this, we note that when the velocity gradients are specified arbitrarily 
as continuous functions of the time ¢, (1) becomes a system of first order 
ordinary differential equations for determining ¢,,(¢). There will exist a 
unique solution taking on arbitrarily prescribed values at any given time fo. 
The velocity gradients can be chosen so that the symmetric tensor d;; takes 
on any desired value at to, so ¢t,; and d;; can be assigned arbitecily at ty. 
Thus (4), must be an identity in ¢,,, Ptbe diem: 
From (2) and (5), we havet 


tiem 


ep " 
¢—~ = D,,;, (6) 

Ct; 
where D;; = tim kmij (7) 
is a tensor invariant of the stress. As such, it is expressible in the form§ 
Di; —_ Ag bj; +Ay tj +Agtixty;, (8) 


where the A’s are scalar invariants of ¢;;. Equation (8) follows from 
(7) whether or not a hypo-elastic potential exists. 


3. We first obtain formal conditions which are necessary and sufficient for 


+ (4), means the fourth equation in (4). 
t In forming derivatives with respect to tijs we replace t;; by its equivalent 4(t;; 


the function differentiated, so that 0/ét;; = 0/ at; ;. 
§ For the case where D,; is a polynomial in the stress, this follows as a special case of a 
result due to Weitzenbéck (9). The A’s are then polynomials in /,, J,, and J;. For more 
general functions, it follows as a special case of the result given in (10, section 40), which is 
derived assuming no continuity. References to other proofs are given in (11, section 6). 


+4) in 
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the existence of a hypo-elastic potential. In N dimensions, the equations 
ov, Ww BY» Ovg ov, : ov Ov, = wal v 
loa oy* ) rye as +i ayr)"8 Been ane 
are satisfied if and only if the vector field v,(y®) is proportional to a gradient, 
v, = ¢ eyb/ey®, (9) being the integrability conditions for the latter system. 
Regarding D;, as a vector in six-dimensional stress space, we obtain from (9) 
(oY Dy_ ot \D nt (Qe Beat (Gt ae) Pam —0, (10) 


Cc ne Cc ti; ot, Cth Ot; at,. 





as necessary and sufficient conditions for the existence of functions ¢ and % 
satisfying (6). To conclude that a hypo-elastic potential exists if and only 
if (10) is satisfied, we must show that the solutions ¢ and ¢ of (6) are scalar 
invariants of ¢;;. Formally, we must show that, for arbitrary ¢ 


(ti) = (tis), W(t) = P(tj;) 
whenever ¢}, is related to t;; by an orthogonal transformation : 
ti = thmAn Amy Api Ans = 45- 
We first consider y. For infinitesimal transformations, we may set 


Aj; = 6;+Q;; = 8;;—Q,, and linearize ¥(t§)—y(t,;) with respect to Q,,, 
obtaining 


ij? 


With) Ht) = © ts) 


= a (tim Quist tim Qi) =2 ob tim Q 


‘mj 


et;; 
= Enj Qiajs 
7 ia é 
where E,, = —-t —— base 
ij ty;  - ote, kj 


Using (8), one can easily show that Z,; = 0 when ¢ is a solution of (6). The 
equations E,; = 0 constitute a system of three independent linear first 
order partial differential equations for y. A routine calculation shows that a 
general solution is an arbitrary function of the functions J, introduced 
earlier, from which it follows that ~ is a scalar invariant of t,;. As a con- 
sequence, @y;/ét;; is a tensor invariant of ¢;;. From (6), 


+= (tata 


from which it is clear that ¢ is also a scalar invariant of t;;. Thus a hypo- 
elastic potential exists if and only if D,; satisfies (10). 
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It should be noted that D,; does not determine ¢ and ¥ uniquely. We 
have , 


oh = 


at;; atti 
~i 
if and only if d= #(iy) , ~ = Fi’), 
. dys 
where F(’) is any continuously differentiable function of ~’ such that 
dF /df’ #~ 0. 


4. Wenow obtain a characterization which seems to have greater intuitive 
appeal. Consider any material for which a hypo-elastic potential exists, so 
that (3) and (5) hold. Suppose a particle P is, at some time t’, subject to 
a non-zero stress ¢,, and that the neighbouring material is deformed in such 
a way that ¢;;d;,; = 0 at P. From (3), either ¢ = 0 or % remains constant 
at P. From (5), 6 4 0 in the neighbourhood of the non-zero stress {;;. 
Furthermore, from (5), not all the derivatives op/ét;; vanish at ¢;;, so the 
equation ¥(t,;) = %(t;) cannot hold for all t;; in the neighbourhood of ¢;;. 
Thus, when (3) holds, there exist, in the neighbourhood of every non-zero 
stress, stresses which cannot be attained by continuously deforming the 
material in such a way that t;;d,; = 0. 

Whenever (1) holds, we have 


Dt;; 
t;; dj; D; ae’ (11) 
where D,; is given by (7). Using (8) and (11), we obtain 
ti d,; = Dyti;. 


If we choose the ¢;; to be functions of the time satisfying the Pfaffian 
equation ‘ 
: Dt, = 9, (12) 


we may assign the quantities w,; -w,; arbitrarily as functions of time 
and use (1) to calculate d;;. Whatever be the choice of w;;, we have, from 
(12), t;;d,; = 0. In (12), the D;; are continuously differentiable functions 
of t,; which, by (2) and (7), all vanish if and only if t;; = 0. In some neigh- 
bourhood of an arbitrarily prescribed non-zero stress, we can choose 
coordinates so that D,, 4 0. Then (12) becomes 


t 


ij 


b+ Fotis t+Fistis + Fos tos + Foo tee + Fas tss = 0, (13) 
where the F’s, given by 
F,;,= D,/Dy when t=j 
= 2D,,/D,, when i ¥), 
are continuously differentiable functions of stress in the region considered. 
In this region, we assume that, in every neighbourhood of each non-zero 
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stress f;;, there exists a stress which cannot be joined to t;; by a solution of 
(13). Then, by a theorem due to Carathéodory (12), there exist, at least 
locally, functions a(t,,, tyo, tyg, tog, tes, tg) aNd B(ty1, tye, ty, tog, tee, gg), Where 
a is continuous and £ continuously differentiable, such that 
op eer es 
l os ol F,=a— (t<j,8 2). 
* tas ye, Sern 
Setting d= Dy, a, Y = Bt, Ht +ter), Htist+tar), H(tes+tsa), too, tys}, we 
obtain ep 
D,; = 4—. 
” *a, 

Our assumption thus implies the existence ot solutions to (6). As noted 
above, a hypo-elastic potential exists when solutions to (6) exist. 

To sum up, we have shown that, for a hypo-elastic potential to exist, tt 
is necessary and sufficient that, in every neighbourhood of each non-zero stress, 
there exist stresses which cannot be attained by continuously deforming the 
material in such a way that t,;d;; = 9. 


5. Results analogous to those derived here can be deduced for the 
situation obtained by replacing t;; d;; by M = (tj; —}ty, 5;;) d;; in (3). These 
may be of interest in plasticity theories of the type discussed by Green (13, 
14), which involve different equations of the type (1) for loading (M > 0) 
and unloading (M < 0). 

An analysis due to Caprioli (15) leads us to conjecture that if the con- 
stitutive equations imply that the work done by the stress in deforming 
any initially unstressed material volume is non-negative, whatever be the 
deformation, then a hypo-elastic potential exists and ¢ can be taken to 
be the mass density. 

The constitutive equations (1) are not appropriate for incompressible 
materials, so our results are not directly applicable to them. It should be 
a simple matter to work out similar results for these, using the constitutive 
equations proposed by Noll (4, section 14). 
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DIFFUSION OF LOAD FROM A BOOM INTO A 
RECTANGULAR SHEET 


By J. B. CALDWELL (Admiralty Offices, Bath) 


[Received 13 December 1956] 


SUMMARY 

A solution is obtained for an elastic plane-stress problem in which load is applied 
to a finite rectangular sheet through a boom riveted to one of its longer edges. The 
other edges of the sheet are stiffened, and resistance to bending in the plane of the 
sheet is provided by an elastic restraint along the loaded edge. Using a stress func- 
tion in Fourier series form, both the efficiency of the sheet under load and the 
maximum shear stress in the structure are obtained as the sums of infinite series from 
which the relative importance of the principal parameters of the problem can be 
deduced. It is shown that the flexibility of the elastic restraint and of the riveted 
connexion have significant effects on the efficiency of the sheet and on the magni- 
tude of the maximum shear stress respectively. 


1. Introduction 

A suip’s superstructure, although it may extend over only part of the 
length of the vessel, must contribute in some measure to the flexural 
strength of the hull, by virtue of its attachment to the upper deck. The 
determination of the additional hull strength in way of the superstructure 
involves essentially the problem of diffusion into a finite rectangular sheet 
of load applied to one of its longer edges. 

Fig. | illustrates the problem to be considered. A load P is applied at each 
end of a ‘boom’ (representing the upper deck) of constant cross-sectional 
area A, to which is attached a rectangular sheet (the vertical sides of the 
superstructure) of length / and breadth 6. The boom and sheet may be of 
different materials, and are connected by a single line of rivets of known 
shear stiffness having constant diameter and spacing. Coordinate axes 
of x and y originate at the centre of the loaded edge of the sheet, and the 
edge y = 6 is reinforced by a ‘stringer’ (the top of the superstructure) of 
constant cross-sectional area a. It is assumed that the shear and flexural 
stiffnesses of boom and stringer are negligible by comparison with those 
of the sheet in its own plane. At the edges x = +4/ transverse displace- 
ments of the sheet in the y-direction are prevented by ‘stiffeners’ (the end 
bulkheads of the superstructure) having infinite extensional stiffness but 
zero flexural stiffness. Lateral displacement of the structure is resisted by 
an ‘elastic foundation’ (representing transverse stiffening of the upper 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958) 
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deck) assumed to be continuous and of constant flexibility along the length 
of the sheet. 
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Fic. 1. Structure and loading 


It is required to find firstly the proportion of the applied load carried by 
the sheet and stringer, and, secondly, the maximum stress in the structure. 
Of the various methods of analysis which have been developed for similar 
diffusion problems, neither the ‘finite-stringer’ method (1, 2) nor the 
‘stringer-sheet’ method (3, 4, 5) is applicable to this case; since the former 
assumes the sheet to be incapable of sustaining direct stress, and the latter 
implies that the sheet is infinitely stiff in the transverse direction and 
therefore cannot be used if the boundary conditions involve transverse 
stresses. Goodey (5) and Mitchell (6) have given more rigorous solutions 
using Fourier integrals to solve the equation governing plane stress distribu- 
tion in semi-infinite sheets. These solutions, and most approximate 
methods, concern the diffusion problems where, unlike the present case, 
both structure and loading are symmetrical. Furthermore, since the 
dimensions of superstructures vary widely according to the type of ship, 
we require a solution that is applicable to sheets of finite dimensions. 

In this paper a Fourier series form for the stress function is obtained 
which satisfies the plane stress equation and the conditions at the edges 
of the sheet. It is shown that the efficiency of sheet and stringer in carrying 
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=n applied load can readily be determined from a series involving explicit forms 
of the principal parameters of the problem. A simple relationship between 
the efficiency series and the magnitude of the peak shear stress in the 
structure is established, from which the conditions for a theoretically 
infinite shear stress are derived. Some general conclusions concerning the 
effects of flexibility of the elastic foundation and the boom-sheet connexion 
are given. 
2. Notation 
a cross-sectional area of stringer. 
A cross-sectional area of boom. 
b breadth of sheet. 
c = cosh ab; c, = cosh ay. 
C constant of integration. 
E, elastic modulus of boom. 
E, elastic modulus of sheet. 
i stress. 
F force in sheet and stringer. 
k flexibility of ‘elastic foundation’. 
l length of sheet. 
n harmonic number (odd values). 
P, force in stringer. 
re force in boom. 
by r applied load. 
<a Q, R, S functions of 8 and A. U = Q/R; V = S/R. 
lar 8 - sinh ab; s, = sinh ay. 
ihe t thickness of sheet. 
_ T function of B,Aando. W = T/R. 
- u,v displacements in x, y directions. 
nd x,y coordinate directions. 
_ x = nr/l. 
_ B = wb = nob/l. 
yu- 1+A+A 
ate Y = am Pe 
=, A determinant involving f, 6, A, A, o, and #. 
he n efficiency. 
ip, E,(u—uz , 12 ine 
6 = (' - , rivet shear flexibility factor. 
b Jax ro 
ed A —a/bt. r = 1422. 
a A = Ap/bt. 


ng - E,/E,. 


~ 
= 
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a Poisson’s ratio. o’ = 1+a; 0” = o—1. 
d stress function. 

m7 = E,kt/b, foundation flexibility factor. 
3. 


Stress function, stresses, and displacements 


The stress function ¢ defining the stress distribution in the sheet must 
satisfy the plane stress equation 


Vid = 0. 
A suitable expression for ¢ is 


— (, nz i 
=> (6, b%eosh("7H) +c, besinh("7#) + 


n=1 


+, by cosh("74) +0, by sinh(“7") oos("7"), (1) 


in which 6 and / are breadth and length of the sheet respectively, C,, , are 


constants to be determined, and m has odd values only. In what follows 
we will consider only the nth term of this series. Putting 


nt . 
a= —; B = ab, coshay = ¢,, sinh ay = 8,, 


T? 
then the stresses and displacements at a point (x,y) in the sheet are 
fe = [Cy Brey +CaBPsy+CyB(28, + aye,)+C, Bl2cy-+ays,) eos ox 
fy = —[O1 Bey +O, B's +Cy aBycy +C, oBys, Joos aux 
fey = [C, B28,-+C; Bey +Cy(Be, +aBys,)-+C,(Bay+aByc,)}sin ox 

E>*[C, 0’ Bbey +C, 0 Bbs, +-C,(2bs, +0’ Bycy) + ee 
+C;,(2be, +0’ Bys,) |sin ax 

v = —E;{C,0'Bbs,+C, 0’ Bbc, +C,(0’Bys,+o"be,) + 
+C,(0'Byc, +0°bs,) eos ax 


u = 


| 
i 





The presence of edge stiffeners of infinite extensional but zero flexural 
stiffness requires that 


f,=v=0, when x= +4). 
These conditions are satisfied by equations (2). 


4. Derivation of the boundary conditions 
The load P applied to the structure can be expressed in the form 


P = > pa cos(“7*) = S <7 sin(dn)eos("7"), (3) 


nodd 





st 
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and we shall deal with the nth term of this series. The first boundary 
condition is that the rate of increase of load P, in the boom plus the load 
transferred to the sheet by shear in the connexion must equal the rate at 
which the applied load increases. That is, if t is the thickness of the sheet 


dP, d 
, + t(fry)e0 —_ ay (Pn cos oat). (4) 


A relationship between the boom load P, and the stresses in the sheet at 
y — 0 is obtained by considering relative displacements of boom and sheet. 
If u, and u are displacements of corresponding points in boom and sheet 
respectively, then due to the shear flexibility of the connexion, uw, will 
differ from u by an amount depending on the shear stress (f,,,),9 at that 
point. We assume that the rivet line may be replaced by an equivalent 
continuous connexion of constant shear flexibility. It is convenient to 
introduce a shear flexibility parameter 6, such that 


whos wag si 


After rearrangement and differentiation equation (5) becomes 


B, 1B — [2 ano “| 
od 2,0 








* dx d. dx 


since @ is independent of x. If the elastic modulus of the boom is E,, the 


load in the boom is 
; dup. 


‘dx’ 


2 = AE, 


and since for the sheet 
E du’ 
8 iz), eae (f-—ofy)20 


the relationship between boom load and stresses in the sheet may be 
written 
dx 


in which » = #,/E,. Substituting equation (6) in equation (4) the first 
boundary condition becomes 


Ap 5 [fe—ahy—0 an 


The second boundary condition concerns transverse displacements of 
the structure under load. If the displacement of a unit length of the ‘elastic 
foundation’ (see Fig. 1) due to a unit load is k, then the transverse dis- 
placement (v),. of the structure is related to the transverse stress (f,),o 


by the equation (v)2.0 = kt fy)20 " 


P, = Au|f.—ofy—0 z2\.. (6) 


+t(fry)z0 = —P, asin az. (7) 
0 


z. 
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since the boom is assumed to have zero shear stiffness. We now introduce 
a dimensionless parameter ¥%, termed the ‘foundation flexibility factor’ 
such that E, kt 


¢= - (9) 





Equation (8) can then be written 
(E,v—wbf,)20 = 9, (10) 
which is the second boundary condition. 
At the edge y = 6 of the sheet the transverse stress in the sheet must be 


zero since the flexural stiffness of the stringer is zero. The third boundary 
condition is therefore (fyen = 0. (11) 


The fourth condition is obtained by considering the equilibrium of an 
element of the stringer. If the direct force in the stringer is p,, then 
dp 
a — "Seven = 0. (12) 


The longitudinal strains in sheet and stringer at their connexion are equal 
(because of equation (11)) and hence if the materials of sheet and stringer 
have the same elastic modulus £,, then 

Pz = a(fz)r» 
where a is the cross-sectional area of the stringer. Using this relationship in 
equation (12), the fourth boundary condition is 


d 
a dz iano" Sew)2p = 0. (13) 


It can be shown that the first and fourth conditions taken together satisfy 
the requirement that the total load in boom, sheet, and stringer at any 
section is equal to the applied load at that section. 


5. Solution of the boundary equations 


Equations (7), (10), (11), and (13) define the conditions to be satisfied 
by the stress function (1). By substituting equations (2) into these 
equations, the following are obtained 


C, Apo’ B® —C,( A6B  +-btB) —C;( ApOB?+-bt)+2C, Aus = P, 
C, 8*—C,0’B—C,0” = 0 
C,c+C,8+C,c+C,8=0 } (14) 


C;Ba-+CyBe+Cy| e+ Ba + ai) +4] 0+B¢(1 + ai) = 0 


in which c = coshab = cosh, 8s = sinhab = sinh®. 








ce 


9) 


0) 
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We now introduce ratios A and A defining the relative ‘effective’ areas of 
boom and stringer to that of the sheet by the pond 
aE = z me Ae i 
Then, putting A’ = 1+ 2A, equations (14) may be written 
Ao'B? = —B(A@B*+1) —(A6B?+-1) en li P,, [bt 


A = 








er ns fl_|e| 
1 tanh B 1 ee 0 | 
BtanhB B 1+A’BtanhB tanhB+ ’B | C 0 | 
or |A||C| = |D}. (15) 
It can be shown that 
» ie  o’d’sech?8—o tanh’s | 9) tanh 
0, = £*\o'x'sechtp— (Se : 2) 
; P.t tanhp  X’ vat 
$f = — r' h? ———_—_—- 
2 = Fpy| YA Baech*B+o (' e tR7 SB si 
C, = cS o’(tanhp— — et F)—vttanh*s-+ 2n8 tanh A)| 
‘ Fat ile 1\ o”sech?8 - 
C= xi le o'(2Atanh B + 4)— y+ oltanh B+) ptanh*p—)| 








in which A = of (A0B?+ 1)Q+AS8]+ (A0B?+1)R+AT (17) 
and Q = tanh?8+ 2A tanh B—X’B*sech?B 


ore ee tanh) 


S = 2(8 tanh B+A’p? tanh®8— p?) 
T = o'*)’B*sech?B+ 20'(3—o)AB tanh B+ 4—o”*tanh?8 

From equations (2), (16), (17), and (18) the stresses and displacements 
in the sheet due to the load P, cos ax can be calculated. By summing the 
effects of each component of the load series the resultant effect of a load P 
acting on the structure may be found. 

The solution indicates that the stress distribution is dependent on the 
dimensionless parameters f, A, A, 6, and ¥, and to a very small extent on 
Poisson’s ratio o. Particular cases of the general problem can easily be 
derived. Putting % = 0 corresponds to an inflexible ‘elastic foundation’ 
which prevents transverse displacements of the boom. This solution 
therefore applies to the symmetrical counterpart of the structure of Fig. 1, 
that is, where the sheet extends equally on each side of the boom. If = 0, 


(18) 
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there is no restraint against lateral displacement. Both limiting cases are at 
present being studied experimentally. A further simplification follows from 
putting @ = 0. This case corresponds to a continuous infinitely stiff con- 
nexion between boom and sheet, a condition more appropriate to a welded 
than a riveted connexion. Some consequences of this simplification are 
discussed in section 8. 


6. Efficiency of sheet and stringer 


If the sheet and stringer were fully effective under the load component 
P,,cosax, then since the longitudinal strain in boom, sheet, and stringer 
would be equal, it follows that the longitudinal stress in sheet and stringer 
would be (P, cos ax)/(Au+bt+a). Hence the total load in sheet and 
stringer at a section distant x from the centre-line of the structure would be, 
in the fully effective condition: 

ee (bt+-a)P,, cos a (1+A)P, cos ax (19) 
' Ap+bt+a 1+A-+A 
in which A and A are the area ratios previously defined. In general, however, 
the load is not distributed uniformly across the width of the sheet. The total 
load F, in sheet and stringer can be obtained by integrating the shear flow 
at the junction of boom and sheet, viz. 


il 
F, —= —t | (fry)z0 dx, (20) 


and from equations (2) and (16) it follows that 
te come (C, B?+-C; B)sin ax 
—P.B , 
= An - | w(tanhe3+ 28 tanh B—A’B*sech?8) +- 


+ 2( REN —tanh*) Sin ax, 


Using the functions Q and R defined by equations (18), 


(fry)z0= — rs —)sin on. (21) 


Then if the efficiency »,, of sheet and stringer under a load component 


P,, cos ax is expressed as the ratio of the load F, to the ‘fully effective’ load 
F’., we find that 


ul 
— (PB/bt\(Q+ Rye | sinax de 
In =F. A(l+A\(1+A+A)P, cosar™ 


~~ 
y 





(22) 
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in which y= ed 
and Mn = were. (23) 


Hence, from equation (17), 
, ~U+l1 
In = (AOB?+ 1)U +AV]+A0p?+1-+-AW 
in which U = Q/R, V = S/R, and W = T/R. 
Finally, by summing the loads in sheet and stringer due to each com- 
ponent P, cosax of the applied load, the efficiency of sheet and stringer 
under the total load P is 





(24) 


= 2% InP, 5 Fi, COB ae (25) 
SP, cos ax 
and using equation (22) for »,, aes ee (3) for the load series 


a <> ts Te gin” cos ax. (26) 


Thus the efficiency of the structure pinowe on the ratio y of cross- 
sectional area of the structure to that of the sheet and stringer and upon 
the terms 7»;, defined by equation (24). Each term in this ‘efficiency series’ 
involves the functions U, V, and W which contain A, 8, and o only. These 
functions are fundamental to the solution of the diffusion problem and have 
therefore been computed. In Fig. 2 values of U, V, and W for0 < 8B <5 
are given, and corresponding expressions which are valid for 8 > 5. 

Some general conclusions can be deduced by inspection of the efficiency 
term 7;,. Since U, V, and W are always positive, as are the flexibility 
factors 6 and x, it can be shown that 7;, decreases with increasing values of 
x; that is, the efficiency of sheet and stringer decreases as the flexibility of 
the ‘elastic foundation’ increases. Similarly, it is evident that efficiency will 
decrease as the rivet shear flexibility factor @ increases, although as will be 
shown in section 8, this will lead to a reduction in the magnitude of the 
maximum shear stress in the sheet. 


7. Maximum efficiency and rate of load diffusion 
The efficiency » of the sheet and stringer is a maximum at the centre of 


length x = 0. Hence 
max = “YS Bac 1), (27) 
7 n 
Values of »,,,, have been computed from equation (27) for structures in 


which both rivet shear flexibility @ and foundation flexibility % are zero 
(Fig. 3), and in which 6 = 0, 4 = o (Fig. 4). In each case three values 
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of the boom/sheet area ratio A, and two values of the stringer/sheet area 
ratio A have been taken. 
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Fic. 3. Maximum efficiency of sheet and stringer. Zero foundation flexibility 


Comparison of Figs. 3 and 4 indicates that if ¢ = 0, i.e. lateral displace- 
ment is prevented, the applied load will ultimately diffuse fully into the 
structure provided the length of the sheet is large compared with its breadth; 
but if ~ = ©, ie. lateral displacement is unrestrained, the maximum 
proportion of the load taken by sheet and stringer is largely dependent on 
the relative areas of boom and sheet. The influence of A and Poisson’s ratio 
oc is generally small; the results for ys = 00 are independent of o. 

Since the maximum efficiency in a particular case is a measure of the 
degree to which the load diffuses from boom to sheet, then the distance from 
the end of the sheet at which diffusion is complete is } times half the value 
of 1/b at which the efficiency reaches its maximum value. Here the effect 
of foundation flexibility is most marked. For example, when 6 = 0, 
A = }, and A = 0, Fig. 3 indicates that if y = 0 diffusion is virtually com- 
plete within a distance of about 3b from the end of the sheet; whereas if 
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ys = o (Fig. 4) diffusion is complete at a distance from the end of the sheet 
approximately equal to its breadth 6. Thus the effect of bending of the 
structure in the ry-plane due to eccentricity of the applied load is to decrease 
the proportion of load carried by sheet and stringer, but to increase the rate 
at which this load is diffused into the sheet and stringer. 
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Fic. 4. Maximum efficiency of sheet and stringer. Infinite foundation flexibility 


If lateral displacements are unrestrained, y% = oo, then as 1/b + oo the 
results of the foregoing theory should approach those of the conventional 
theory for tie bars under eccentric end tension. The latter indicates a non- 
linear response to applied load, arising from the term sech,/(P/?/4B), in 
which B is the flexural stiffness of the structure referred to bending in the 
xy-plane. The present theory, on the other hand, is linear with respect to P. 
It can be shown that if //b = oo, the maximum efficiencies predicted by the 
two theories are identical if sech,/(P/?/4B) = 1, and since for the type of 
structure under consideration ,/( P/?/4B) is unlikely to exceed 0-3, the dis- 
crepancy between the theories will always be less than 5 per cent., so that 
the linearity of the present solution does not invalidate the results obtained. 
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8. Maximum shear stress in the sheet 


It follows from equations (3) and (21) that the shear stress at the junction 
of boom and sheet is 


Fenda = ~ > ("F) 5) oS sinew /2)sin(ax), 


and therefore, using equation (23), 
(fry)2,0 — wd a Nn sin(n7r/2)sin(az). 


In particular when x = }l, since sin*4n7 = 1 for odd values of n, the 
maximum shear stress in the structure is 


4P , 
(fry)iro - Uahou seal atin >, Tn: (28) 


This relationship enables the peak shear stress to be calculated from the 
efficiency series—equation (24); and furthermore the condition for a 
theoretically infinite shear stress is obtained by examining the convergence 
of the 7, series. Using equations (17) and (18), it can be shown that 
“f w/A 

HOB (GE 2OR 
Then if @ is non-zero, since yw, 6, and £ are all positive, each term in the 7}, 
series is less than the corresponding term of the series 

,  wPlA | 
"mn = Sop? — NOB 


But since this latter series is known to be convergent (for n odd 


(1), )p—a 


> 1/n® = 4n*) then the series for 7),, equation (24), is convergent. Therefore 
for non-zero values of the rivet shear flexibility factor 0, the maximum 
shear stress is finite. 

If, on the other hand, the connexion is infinitely stiff in shear, then 
6 = 0, and the series for 7}, leads to 


‘ 1 
(Mn )p—-a —- Ap’ 


which represents a non-convergent series. In this case the peak shear 
stress is theoretically infinite. 

The effect of rivet shear flexibility on the maximum shear stress is 
illustrated in Fig. 5, where for y = 0, A = 1, and l/b = a, the ratio of 
maximum shear stress to maximum direct stress f, in the boom clear of the 
sheet is given for three values of A and 0-01 < @ < 100. Tests on single 
riveted lap joints by Flint (7) indicate that in ship structural connexions 
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0-01 < @< 0-1. Within this range the maximum shear stress is unlikely 
to be excessive unless the boom area ratio A is large, or the modular ratio 
p = E,/E, is small. 
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Fic. 5. Maximum shear stress in sheet 


9. Conclusions 

A Fourier series form for the stress function can be used in the two- 
dimensional elastic problem of diffusion into a finite rectangular sheet of 
load applied to a boom riveted to one of its longer edges. It is shown that 
lateral resistance to bending of the structure in its own plane has a significant 
effect on the rate at which load is diffused into the structure, and upon the 
efficiency of the sheet in carrying direct load. A relationship between this 
efficiency and the maximum shear stress in the sheet has been established, 
from which it has been shown that the peak shear stress varies with the 
shear flexibility of the connexion between boom and sheet, and is infinite 
when the rivet flexibility is zero. The effect of variations in the geometry 
and elastic properties of the structure can be determined without difficulty. 
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SUMMARY 

Consideration is given to highly elastic materials which may deform subject to 
a three-dimensional system of constraints, such as would be produced by three sets 
of ideally thin, flexible inextensible cords whose paths are defined by the coordinate 
curves of a general curvilinear reference frame related to points in the undeformed 
body. When the cords of each set lie initially in parallel straight lines, the equations 
governing the deformatior may be solved in terms of arbitrary functions. The 
nature of the solution for the general case is also examined. 


1. Introduction 

In the theory of finite elastic deformations, some consideration has recently 
been given to materials which may deform subject to certain types of 
geometrical constraint (1, 2, 3, 4). Examples of this feature are provided 
by incompressible materials, where there are no volume changes during 
deformation, and by bodies reinforced by systems of ideally thin, flexible, 
inextensible cords, which prohibit extension in certain directions in each 
element containing them. Each such constraint thus implies a geometrical 
relationship between the variables defining the deformation and introduces 
an arbitrary parameter into the stress system. 

A case of special interest arises when a uniform plane thin sheet of elastic 
material is reinforced with a network consisting of two sets of thin, inexten- 
sible cords lying in its middle plane. This problem has been considered by 
the writer (2), and the very similar problem, of a network of cords, without 
elastic material, by Rivlin (8). In both cases a pair of differential equations 
is obtained for the determination of the deformation, and if the cords lie 
initially in parallel straight lines, these equations may be solved in terms 
of arbitrary functions. It may be shown quite generally that the system 
of differential equations for the determination of the stress and displacement 
functions is hyperbolic in character. 

In the present paper the analogous three-dimensional problem is examined 


+ Now at the University of Nottingham. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958) 
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briefly, it being assumed that any deformation is subject to the restriction 
that there shall be no extension along the coordinate lines of a curvilinear 
reference frame related to points in the undeformed body. The equations 
for the determination of the stress and displacement functions are again seen 
to be hyperbolic, the characteristic surfaces being those defining the system 
of constraints; when these surfaces are planes, the general solution can be 
given in terms of arbitrary functions. 


2. Notation and formulae 

The notation employed is that originated by Green and Zerna (5, 6) and 
extended by the writer (2). The points of an unstressed and unstrained body 
at rest at time ¢ = 0 are defined by a system of rectangular cartesian 
coordinates x‘, or by a general curvilinear system of coordinates 6‘. The 
latter reference frame moves with the material as it is deformed and forms 
a curvilinear system in the strained body at time t. Points of the deformed 
body may also be defined by a set of rectangular cartesian coordinates y’, 
and in the present paper we shall take the 2z‘-axes and y-axes to coincide. 
With the coordinate system @ in the unstrained body we may associate 
covariant and contravariant metric tensors g;;, g"’, the corresponding quan- 
tities for this coordinate system in the strained body at time t being G,;, G. 
The covariant strain tensor y,;, referred to the coordinates 6‘, may then 


be defined by 
e defined by Yu = UGy—Gy)- — - (2.1) 


For a second system of curvilinear coordinates 6’, metric and strain tensors 
9i;, 9°, G,;, @", and 7,; may be similarly defined and may be derived from 
the corresponding unbarred quantities in the usual manner by tensor 
transformations. 

We consider here elastic materials subject to constraints which may be 
described by means of functional relationships 


SnlFg) = 0 (m = 1, 2,...,n, n < 6) (2.2) 


between the strain components 7,;, the restriction upon the value of n 
resulting from the consideration that more than five such constraints 
would permit only rigid-body motions. For such materials it has been 
shown by the writer that the stress tensor 7) referred to the convected 
coordinates @' in the deformed body, may be written in the form 


i — ti MS 2 Fm (28 08) , 20" a6) : 
vn re +4 2, Ym cae age t abe 6) (n < 6), (2.3) 


where g,, are arbitrary scalar functions of 4’, and 7// is a function depending 
upon the constraints only to the extent to which the relations (2.2) introduce 
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simplifications in the form of the strain-energy function W. The exact form 
of 7’ is unimportant for the purpose of the present paper; the general 


expression is (2, 6) Pins 1 fig aw ow 


‘. ee j _ - - 
where J = (9431, G= |G;|- 


In the absence of body forces the equations of equilibrium may be 
written ril|, = 0, (2.4) 
where the double line signifies covariant differentiation with respect to 
the deformed body, that is, with respect to 6‘ and the metric tensor com- 
ponents G,;, G". 


3. The constraint conditions 


The two-dimensional problem in which a plane uniform elastic sheet is 
reinforced in its middle plane by means of two sets of thin, flexible, inexten- 
sible cords following the coordinate lines of a plane curvilinear reference 
frame, and is subjected to a plane deformation has been examined by 
Adkins (2). The analogous problem for a plane network of initially 
straight cords, without elastic material, has been considered by Rivlin (8). 
For the case of an elastic sheet, it is assumed that the cords restrict the 
deformation so that there is no extension along any member of the families 
of curves in which they lie. 

A natural extension of this problem to three dimensions is obtained by 
postulating a material which is such that in any deformation there can be 
no extension along any of the coordinate curves of the reference frame @. 
If corresponding elements of length in the body before and after deforma- 
tion are denoted by ds, and ds respectively, we have 


dst = g,,d0'd0) = g,,d0‘d6i, 
ds? = G,;,d0'd0' = G,, d0‘dd', 


and since ds, = ds along the curves 6‘ = const, 6’ = const (i + j), these 
expressions lead to the constraint conditions 


Gi = G, or Yq =O (i not summed), (3.1) 
which may also be written 


sas as i 
Re ies Se es ee wna i not summed 3.2 
«= op op = a= sp gp not summed). i 


If the coordinate surfaces 6‘ = constant in the undeformed body are 
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planes, so that the components of the metric tensor g,;; are constants, 
equations (3.2) may evidently be integrated to yield 

y" = f,(O) +fo(F) +fo(®) 

y® = 91(2)+-92(8)+-99(%) }, (3.3) 

y= h, (8) +-ha(8*) +hg(®) 
where f,(), g,(0"), and h,(8") are functions whose first derivatives, denoted 
by primes, satisfy the relations 

SP+GP +h? = Gr (3.4) 


but which are otherwise arbitrary. 


4. Equations of equilibrium; stress strain relations 
For the purpose of this paper, it is convenient to introduce a symmetric 
covariant tensor ¢,; which is such that 
ii x thm lan, (4.1) 
where e'/* takes the value 1/VG or —1/VG@ according as i, j, k is an even 
or odd permutation of the numbers 1, 2, 3, and is zero if any two indices 
are equal. It is easily verified that when r’ is expressed in this manner, 
the equations of equilibrium (2.4) are identically satisfied. Moreover, this 
remains true even if the off-diagonal components ¢,; (i  j) are zero. We 
therefore choose ¢,; so that the components ¢,; for i # j are zero in the 
reference frame 6! defining the system of constraints. From (4.1) we then 


— FH BMI ny, (4.2) 
where the double line now denotes covariant differentiation with respect 
to 6‘ and the metric tensor components G,;, G/. The components ¢,, 
correspond to the stress functions of the classical theory of elasticity, 
given, for example, by Love (7). 

The constraint conditions (3.1), when combined with (2.3), yield for + 


the formulae 3 5 ant 
io * 06* 60) 
1} — 7 > we 
ri = 4 pa Qe abe ole’ (4.3) 
or, if the coordinate systems 6‘ and 6 are allowed to coincide, 
3 
Fi = FH > q, 8. 84, (4.4) 
r=1 
5! being the Kronecker delta. Combining (4.4) and (4.2) then gives 
Fit — Fiitq, = Etkmgikmg, ||, (i not summed), (4.5) 
Fi = FY = ema S linn (bt FJ) (4.6) 


If the deformation has been determined from the constraint conditions 
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(3.2), the components 7 may be regarded as known functions of 6‘ and 
(4.6) furnishes three equations for the determination of the functions ¢,.,.. 
The parameters q; are then given by (4.5); if the constraints arise from 
systems of thin flexible inextensible cords, these parameters may be 
expressed in terms of the tensions in the individual cords and the spacing 
of the adjacent cords of each system. 
When the coordinate surfaces 9 in the undeformed body are planes, we 
see, from (3.3), that the metric tensor components G,, are given by 
Gi = Sift GG hi hy, (4.7) 
so that the derivatives @G,;/20" are zero if i, j, and k are all different; it 
follows that the Christoffel symbols I‘, formed with the barred quantities 
G,;, G® vanish if r ~ s. Equations (4.6) then yield 
an 1 we 1 1 eae 
G eee’ = * G ebee8 G eb\ab? 
Since the components 7!’ are, for a given type of material, known functions 
of f,, g,, h, and their derivatives, we may obtain from (4.8) by integration 


$1, = X03 +0, (8, 8)+ 0,(8, 8) 
$o2 = X13 + O,(F, &)+ 0, (82, ) }, (4.9) 
$3 _ Xie tH (8, B) + (8, 6) 
a gf 
where Xre = -I db | Gz de (4.10) 


are known functions of #', and @,, ®,, and V, are arbitrary functions of 
the arguments indicated. 

The remaining components 7” of the stress tensor are obtained by com- 
bining (4.9) with (4.5) and evaluating the covariant derivatives. Since 


it he trl 
VG=|9 o 9 |, (4.11) 
in g& Ri 
it follows that hae oa es 
GW — ] {onG ovG * é VG evG onG OvG| 
G\ éf' of; og, @ eg; en, oh; u 
and the non-zero components of the Christoffel symbols I'j, take the forms 
1 eG 
2G 6’ 
analogous formulae holding for Tj, and T'i,. Combination of these results 
with (4.7), (4.11), and (4.12) yields 


, onG gf ONG , evG - 
i acage ol node G. 4.14 
n= (l a ag ae) iia 


(4.12) 


Mi, = G86,,,+G4°6,,,, n= (4.13) 
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From (4.5) we may now obtain 


i s ‘aso oe) 1 ‘so oR) 


VG eB\VG 0)" VG ab'\vG ab 
on a m=, @ 1 
Tis cg G)— Mca’) (Po Spt + Te) 


with corresponding expressions for 7? and 7%. 


5. General nature of the equations 

In the two-dimensional problem examined in the earlier work (2, 8) the 
differential equations which enter into the analysis are found to be hyper- 
bolic, and it is then possible to prove that specified distributions of stress 
and displacement along a plane arc are sufficient to determine the deforma- 
tion and the stress system within a region bounded by the given arc and 
the characteristic curves through its ends. The general similarity of the 
basic assumptions and the resulting equations in the present instance 
suggests that analogous results may be obtained for the three-dimensional 
problem. A full analysis for all possible cases would evidently prove much 
more lengthy than the corresponding treatment for two dimensions. We 
shall therefore concentrate attention on proving the following proposition: 

If the displacement components, and hence the quantities y‘, are given 
over any surface S, no part of which coincides with a characteristic surface, 
then the deformation is completely determined within a region V mapped 
out by S and the characteristic lines through its boundary C. 

We may, without loss of generality, choose S to be the surface 6* = con- 
stant, the lines @* = constant (a = 1, 2) then forming a curvilinear system 
on this boundary. For brevity, we use the notation 
= S, ge = =, 
™" oo 7 a 06 
it being understood that in the subsequent work Greek indices take the 
values 1, 2. Since the components y‘ are given over the boundary S, their 
derivatives y', may be regarded as known quantities over this surface, 
whilst the quantities 6‘, are also known from geometrical considerations... 
It is then possible to show that the constraint conditions (3.2) enable the 
derivatives 7‘; to be calculated over S. 


5: — 


(5.1) 


From the relations y?, = 79.4%, 


which hold over S, we may obtain 
Akigh, = Bri, (5.2) 
where Bri = y2 B,—y2B, (5.3) 
At = 5 8,—05,8, (A* = 0))’ 
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and A* transforms as an antisymmetric contravariant tensor with respect 
to transformations of coordinates 6‘. Equations (5.2) yield 


Bri ASigh — Agr, Br? 482g — ANGn 
from which, by squaring, summing over p, and making use of the constraint 


3 3 
conditions (3.2) to eliminate ¥ (y4)? and > (y%)?, we obtain 
p=1 


pl 
3 3 
>, (BP)? — & CPt (Aaa = (A") Gow 
= 
3 
2, (Br2)? — 3 ps + (A**)?g5, = (A™)*Gu1, 
p= 
3 3 
or >> Cn 9's = Xv 2. Cro 9s = Xa (5.4) 
. = 


Xi» Xz and C,, being given by 


3 
ps es bY (.BP1)?+-(A81)?9,,—(A™)*Gos 


os . (5.5) 
mo z, ( BP?)? + (A*)*g5,—(A™)*9, 
p= 
and Ciq = 2BP* A (a not summed). (5.6) 
Writing Ci2X%1— Cu Xe = ‘th, C22 X1— Cun X2 = Vo, (5.7) 
Aye = C41 Cogp—C gC, = 44°1A9*( BU B22 — B!2 B21) 
Aas = Cy_ Ca —Cy, Cog = 44*1A22( Bi? B31 — Bi B82) }, (5.8) 


Ags = Cyq C31 —Cy, Cag = 4A4°1A5*( B22 B31 — B2) B32) 
we may obtain from (5.4) 
I's = (F.—Ags 75)/ Are, Fs = —(H—Ais F%s)/ Are, (5.9) 
and when these expressions are introduced into the third of the constraint 
conditions (3.2), we have 


Jag = {(Azg+ As + A$s)(9%3)?— 2(Ays ‘Vy + Aggy ‘Y2)9°, + ‘V2 + ‘¥3}/A2, 


Solving this equation we obtain 


7°, = {M+(M?2—LN)/L, (5.10) 
where L = Aj,+A?,+A3; 
M = Aj, ¥,+Aes "i (5.11) 


Ty on OE £2 5 2 
N = ¥?+" 3—Iss Ais 
are quantities whose values are known over S. 


The sign occurring in the expression (5.10) is determined from the con- 
sideration that the functional determinant 


vq — yy ¥*) 
——— @(B1, 8, 3) 
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must be positive. For, when the relations (5.2) are used to eliminate the 
derivatives 7%, 7% from the expression for VG, and the resulting formula 
simplified by means of (5.8), we obtain 
VG = (Avs J!3— Ags F3—Ay2 9%s)/ (4A? AA A®), 
Further reduction by means of (5.9) and (5.11) yields the result 
VG = —(Li*,—M)/(4A12A*A51432A,,), 
or, by virtue of (5.10), 

VG = F(M*— LN)t/(4A%A451482A,,), (5.12) 
the upper and lower signs in (5.10) and (5.12) corresponding. All of the 
quantities entering into this expression are determined uniquely by the 
boundary conditions over the surface S and the geometrical parameters 


of the undeformed body, and only one of the signs occurring in (5.12), and 
hence in (5.10), is therefore possible if the expression for VG is to be positive. 





Fic. 1 


When 7*,; has been determined, unique values for the remaining first 
derivatives of y' follow from the linear relations (5.9) and (5.2). Higher 
order derivatives of y‘ with respect to the coordinates 6* may be obtained 
by differentiation of the expressions for those of the first order, the analysis 
following the lines indicated for the two-dimensional problem (2) By 
Taylor’s theorem, an expansion for the coordinates y* may then be obtained, 
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valid within a finite region V bounded by S, and satisfying the constraint 
conditions within V and the boundary conditions on S. By this means, 
the values of y' may be determined over a surface S’, close to but not 
coincident with S; the region throughout which the solution is known may 
then be extended by a repetition of the procedure, starting from S’. 

The region throughout which a solution can be obtained by this method 
may be mapped out as follows. Let characteristic surfaces 6‘ = constant 
intersect the surface S in the curves PQ, QR, RT as shown in the accom- 
panying diagram, the points P, Q, R, and T lying on the boundary curve C 
of S. If PQ, RT intersect in O, the curvilinear tetrahedron bounded by 
the portion OQR of S and the characteristic surfaces through OQ, QR, 
and RO represents a region V’ throughout which the values of y‘, and 
hence the deformation, may be determined by the foregoing process. In 
crossing the surface 6! — constant from a point just inside this region, the 
quantities y', and hence their tangential derivatives 7',, 7',, must be con- 
tinuous, but this consideration does not apply to the normal derivatives, 
and hence to 7‘, which may be assigned any values consistent with the 
first of the constraint conditions (3.2). Similar considerations apply for 
each of the boundary surfaces 6‘ = constant. To determine the values of 
y' outside the region V’, it is therefore necessary to make use of boundary 
conditions on S outside the curvilinear triangle OQR. Regions adjacent 
to V throughout which the deformation may be determined may be mapped 
out by drawing further characteristic surfaces adjacent to the boundaries 
of V’. The dotted lines QO’, RO’ represent the intersection of two such 
surfaces with S. By continuing this procedure until the entire area S 
bounded by C has been covered, the region V may be defined throughout 
which the deformation may be determined by the boundary conditions 
on S. It follows from this analysis that the constraint equations (3.2) are 
hyperbolic, with the surfaces defined by the system of constraints forming 
characteristic surfaces; examination of the equations (4.6) for ¢,, shows 
that these are similar in character. 

If the characteristic surfaces 6! are initially planes, the constraint con- 
ditions have a solution of the form (3.3). In this solution, the functions 
t,, Ir, 2,, are determined apart from arbitrary additive constants, by the 
boundary conditions, and consideration of the second and higher order 
derivatives of the coordinates y' becomes unnecessary. 


This work forms part of a programme of research undertaken by the board 
of the British Rubber Producers’ Research Association. 
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SUMMARY 


In the present paper, a unified theory of internal ballistics of orthodox, leaking, 
recoilless, high-low pressure, recoilless high-low pressure guns, and of solid-fuel 
rockets has been presented. The internal ballistics of the most general weapon— 
R.C.L. H./L. gun—have been studied in detail. 


1. Introduction 

THE main problems of internal ballistics have been almost completely 
solved for the orthodox gun (Corner (1), H.M.S.O. (2)). For the newer and 
still unorthodox guns like the R.C.L. and the H./L. guns, solutions have 
also been given (Corner (1, 3, 4), H.M.S.O. (2)). A theory of the internal 
ballistics of solid-fuel rockets also exists (Crawford (5), Trubridge (6), 
Wimpress (7)). But most of these theories have so far been developed 
independently. A unified theory for all these is still awaited, though 
indication of its necessity was given as early as 1950 by Corner (1) when 
he predicted the union of gun and rocket theory covering in one style of . 
treatment the many new weapons. It is obvious that such a theory cannot 
be a simple one; nor for that matter, it is likely to throw additional light on 
the understanding of the extreme cases—the orthodox gun and the solid- 
fuel rocket—which are already quite well understood. The attempt is 
worth while from the point of view of development of ‘half and half’ 
weapons, i.e. weapons which have partly the characteristics of guns and 
partly those of rockets. 

In the present paper we have made an attempt to develop such a unified 
theory of internal ballistics. The most general weapon considered here is 
the recoilless high-low pressure gun, and we first write down the system of 
equations for this gun. From this we deduce as particular cases, the system 
of equations for R.C.L. gun, for H./L. gun, for orthodox gun, and for a 
solid-fuel rocket. ? 

Before the charge is burnt, the isothermal model for R.C.L. H./L. gun has 
been studied, but since after burning, the gases cool off considerably due to 
expansion, the variation of temperature has been taken into account for 
this case. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958] 
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2. Fundamental equations for a recoilless high-low pressure gun 

The figure below gives the basic diagram for such a gun. 

Let F, C, D, B, ny, & denote respectively the force-constant, mass, 
ballistic size, rate-of- burning constant, co-volume per unit mass, and density 
of the charge. Let K,, K, be the capacities of the two chambers when the 
shot is in its initial position, and let A be the area of cross-section of the bore. 


First chamber Second chamber 
pe ee Te pes! =>) 
3S Ty Ny LS, CS Spon Aaa es 
AB notin Te ner, 
Exit nozzle Middle nozzle 
Fic. 1. A R.C.L. H./L. gun 














At time ¢, let x, z, f denote respectively the shot-travel, the fraction of the 
charge mass burnt, the fraction of the web size remaining, and let p,, 7, p;, 
ON,; Pe, Te, pe, CN, denote respectively the pressures, temperatures, gas 
densities, and gas masses in the two chambers. Then we have the following 
basic equations for the R.C.L. H./L. gun: 


Equation of state for the first chamber 


r C 1— Yar y 
p,| x, CN, | ~ CN, RT, (1) 
Equation of state for the second chamber 
p,| K, + Ax—CN, n| = CN, RT,. (2) 
Equation of continuity for the first chamber 
oF — COM, YSiPi, PSs Pr (3) 


dt ' (RT,)* ' (RT,) 


Here S,, S, are the throat areas of the two nozzles and 


9 (y+D/Ay-D 
= yt{ 4 
lata (Fi) (*) 
Also we have assumed that 
2 \viy-) 
Pe < (=| , : (5) 
Pr yt+l 


If p,/p, is greater than this critical value, y has to be replaced by 


Nee Or ws 


In practice, however, (5) will be satisfied and we shall assume in the rest 
of the paper, unless stated otherwise explicitly, that y is constant. We 
then have 
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The equation of continuity for the second chamber 


dN, ; WS, Py 
Cua > (RT) °) 


so that (3) can also be written as 


dz dN, dN, , S\p 
OF = C=1+0° ve 7 
dt dt + dt * (RT)! (7) 
The equation of motion of the shot 
ae », dx , dx 
(1-050 +30) 73 — Wea =< Apyg. (8) 
The rate-of-burning equation 
d 
D :/ = — Blp,). (9) 
The form-function z= ¢d(f). (10) 


Energy equation for the first chamber 
(i) In time dt, a mass C dz with internal energy C dzC,,7, enters the gas. 
Let 57, be the change in temperature of the first chamber due to this 
additional gas, then 
C dzC,T = CN,C, 8T,+C dzC,T, 
. dz 
or 67, = —(%—T,). 
N, 
(ii) In time dt, a mass C(dz—dN,) escapes through the nozzles. Since 
the expansion is adiabatic, we have if A7;, is the change in temperature due 


to this effect AT, a y-l Ap, _ yr! dN, —dz 


T, \l—mp pi = 1—m, 





Henee, 7 = (y— n+)", 
where eo . ae 
1—mp; 
The change d7, in temperature, due to both the factors, is 
aT, = 8T,+AT, = x to T) + OK 1-+e)(dN,—dz). 

d i. d! N, 

~(N,7,) = T%—— —*}. 
Hence 7 TT) = T 7 [y+ wae = 


Using (3), 


d a dz We 
SMT) = %S—Ly y+ Dd [BAP cay (RT). (14) 





6) 


8) 


9) 


0) 


ce 
ue 
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If we neglect the co-volume term, the energy equation for the first 
chamber becomes 


(M7) = 7, —Ps (ary EPpy, (12) 
€ 

This derivation of the energy equation for a R.C.L. H./L. gun is based 
on Corner’s (3) method of deducing the energy equation for a R.C.L. gun. 
Alternative derivations of the same equations are obtained by using the 
principles used in Kothari (8) or Kapur (9). 

Energy equation for the second chamber 

In time dt 


(i) the addition of internal energy to the second chamber due to the 
energy brought by the mass C dN, entering from the first chamber is 


can( . i) an; 


(ii) this gas, when it forces its way from the first chamber does work on 
the middle nozzle and this energy is also added to the energy of the second 


chamber. This is C dN, RT(1+<); 
(iii) increase in internal energy of the second chamber is 
“s j UM T); 
y- 
(iv) work done on the shot is 
Ap, dz. 
Hence the energy equation for the second chamber is 
= nT + CaN, RT,(1+e¢) = a d(N, T,)+ Ap, dz, 
y—1 
or taking heat losses into account, 
i _A dN, 
T(t) = —Sy—-nS+f+—na+en 2. (1s) 


The equations (1), (2), bis (6), (8), a (10), and (13) constitute the most 
general system of equations, from which the system of equations of all 
weapons will be deduced as particular cases. This we proceed to do in the 
next section. 


3. Particular cases 
1. R.C.L. H./L. gun. Isothermal model. Before charge is burnt 
In this model we assume the temperatures of both the chambers to be 
constant and to be the same. These two assumptions replace equations (11) 
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and (13). In fact if y and 7 are put equal to unity, (11) reduces to (3) and (13) 
becomes an identity. Denoting R7, and RT, by A, equations (1), (2), (3), 
and (6) become 


ps] Ki a — CN, "| = CN,A, (14) 
p,[ K,+Ax—CN, n] = CN,A, (15) 
dz dN, WS,p, . Sep 

Cy = ee 7 ana) et | 21 
a>" ats * 2s” (16) 


IN, wWS,p 
c= .. Pas. 17 
dt ! (17) 
The seven equations (8), (9), (10), (14), (15), (16), and (17) constitute the 
system of equations for this model. The integration of this system will be 
discussed in section 4. 


3.2. R.C.L. H./L. gun. Non-isothermal model. After charge is burnt 
After all the charge is burnt z = 1, f = 0 and equations (1), (3), and (11) 


are modified to p;[K,— CN, ] = CN, RT, (18) 


,4N, . WS, p, , oS, p 
"hye 1 1 Cas 2F1 ae 0, 19 

dt (RT,)' © (RT,)! a 
d 


a bS. Sy p . 
ae T,)+ly+(y— Vel CREAT + PRT) = 0. (20) 


The seven equations (2), (6), (8), (13), (18), (19), and (20) constitute the 
system of equations for this model. The integration of this model for the 
case 7» = 0, « = 0 will be discussed in section 5. 


3.3. H./L. gun. Non-isothermal model. Before charge is burnt 


first chamber 


QW AR NN rd chamber 
Ky Ko \ 


a ee 
— ee ee 
Fic. 2. H./L. gun 


























This case is obtained by putting in the general system of equations 
8, = 0, N, = z—N,, 


so that (3) and (6) become identical and equations (2), (3), (11), and (13) 
are modified to 


p.| K+ Ax—C(z—N,)y] = C(z—M,) RT,, (21) 
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dz _ GAM, MPs 
Ca Oat (RT (22) 
i d S. 
SMT) = To — Lyte] EP (RT, (23) 
d » Ap,,- d. d 
“le-M Ih] = —GRO-NS +y+—DINGe—M)- (24) 


The eight equations (1), (8), (9), (10), (21), (22), (23), (24) constitute the 
system of equations for the model. Its integration for a tubular charge 
with a linear law of burning has been discussed by Kapur (10). 


3.4. H./L. gun. Semi-non-isothermal model. Before charge is burnt 

In this model, the temperature 7; in the first chamber is assumed con- 
stant, while temperature 7, is assumed to vary. The first assumption 
replaces (23). Denoting the constant RT, by A (24) is replaced by 


ds N\(RT,—A)) — —AP2(5—1)" 
al? N,)( RT, A)] = CRY Na. (25) 


The seven equations (8), (9), (10), (14), (21), (22), and (25) constitute 
the system of equations for the model. From (8) and (25) 


diy. ase 
al (z—N,)(RT,—A)] = —(¥ IW aa 


(26) 

This is similar to the form of the energy equation derived by Aggarwal 
(11). He also discusses the integrations of the system for the linear law 
of burning. 


3.5. H./L. gun. Isothermal model. Before charge is burnt 


In the model 7, and 7; are assumed constant and to be the same. These 
two assumptions replace (23) and (24). Denoting RT, and RT, by A, (21) is 
modified to 

p,| K,+ Ax—C(z—N,)n] = C(z—N,)A. (27) 

The six equations (8), (9), (10), (14), (22), and (27) constitute the system 
of equations for the model. These are identical with the equations first 
deduced by Corner (1, 4). Corner discussed the integration for this system 
for the linear law for a tubular charge. The integration for the more general 
form-function (10) has been studied by Aggarwal (12, 13) and Kapur (14). 
Its integration for the general linear law has also been discussed by 
Kapur (14). 
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3.6. H./L. gun. Non-isothermal model. After charge is burnt 
In this case z = 1, f = 0, and (21), (22), (23), and (24) are modified to 


p[K,+Az—C(1—N,)y] = C(1—N,) RT, (28) 
aN, _ _ 9S: (29) 
dt C(RT,)’ 

st a ar a Ps (RT) (30) 

d ™ A a. 

qt Mt] = aR! nS —[yt+l(y— I)e]T, ~ = (31) 


The six equations (8), (18), (28), (29), (30), and (31) constitute the system 
of equations for this model. This model was first discussed by Corner (4), 
but he overlooked the energy entering from the first chamber into the 
second. Consequently his energy equation corresponding to (31) did not 
contain the last term in the right-hand side of (31). The corrected equation 
was used in the discussion of the same problem in Kapur (10). 


3.7. Recoilless gun 


» Ss 
b,T.p,N  +-----#---ld)-----{Shot > 


~ J 


Fig. 3. R.C.L. gun 














In this case 
S=8, n=p=p, T=T1=T7T, p=pr=p, M+N,=N. 
We note that since p, = p, from (4a), the rate of mass flow from one 


chamber to the other vanishes and the term corresponding to this vanishes 
from (3). From (1) and (2) we get 


p|K+4z—“U—)_ own] = ONRT, (32) 
from (7), oS - at: — ‘gi (33) 
and from (11) and (13), 
d 7@ a se »_ Ap(y—V) I) dx . 
aAT) = Te eis saF f ~ 


The six equations (8), (9), (10), (32), (33), and (34) constitute the system of 
equations for a R.C.L. gun. These are the same as the equations for a 
leakirg gun obtained by Corner (3) except that instead of (32) Corner gives 


o| K4 — = 











|= ONRT(1+ (33 a) 
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The difference with this equation, as well as his definition of W,, is due to 
his taking into account Lagrange’s relation for the variations of pressure 
inside the chamber. The integration of this system of equations has been 
discussed by Corner (3). 


3.8. Orthodox gun 








f i 
: Atos ie 
ys seeps 








Fic. 4. Orthodox gun 
In this case 
S=9 yn=pha=p TY=hZ=T, »=p=p, += N=z 
Also since p, = p, from (4a), the rate of mass flow from first chamber 
into the second vanishes and the corresponding term in (3) disappears. 
From the equations (1), (3), (11), and (13), we get 
] 


p| K+ Ax— + cx 


; ")| — CzRT, (35) 


d _md ,._,,Apdz 
ae) = FoR a (38) 
From (8), (35), and (36) 
dj [, C a me x, dx d*x 
. + Ax — _ = —_—(y—1 oe prone 
nA? | wait ox; »)]] CP a 9 "a ap’ 
which on integration gives 
CF: = p| Ate +)—C-(n— 5) ]-440G—DMe* (37) 
il 
where Al=K —5 


and v denotes the velocity of the projectile. The four equations (8), (9), _ 
(10), and (37) give the standard equations for internal ballistics for an 
orthodox gun. 


3.9. Solid-fuel rocket 


Warhead - - m( = 








Rocket motor Nozzle 














Fic. 5. Solid-fuel rocket 


_ S,= 8, S, = 0, Pi = P2 = P, T, = T, = T, 


Aaa N,+N, = N, x= 0. 
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Except for x = 0 the conditions are the same as those for a R.C.L. gun. 
Thus the equations for a solid-fuel rocket can be deduced from those of 
section 3.7 by putting z = 0 in (32) and (34) so that 


p[K—< ; =) _ CN» = CNRT, (38) 
| dz ybSp 

NT) = 2 39 

7“ (NT) = 1 OR (RT) (39) 


The five equations (9), (10), (33), (38), and (39) constitute the system of 
equations for a solid-fuel rocket. The connexion between the theory of 
recoilless guns and solid-fuel rockets was first noticed in Kapur (15). The 
system was also integrated, both before and after the charge is burnt, with 
special reference to obtaining the equation of the pressure-time curve. 


4. R.C.L. H./L. gun. Isothermal model before charge is burnt 

Equations (8), (9), (10), (14), (15), (16), and (17) constitute the system of 
equations for the model. We consider their integration for the usual linear 
law of burning, so that (9) becomes 


df 
Do = — Bp. (9a) 


4.1. Integration of the equations. Maximum pressure in the first chamber 


From (9a) and (16), 
dz_ dN, _y df y df 


a a a ae ©) 
where Y= ave Y= ane (41) 
are the two dimensionless leakage parameters. Integrating, remembering 
that initially z = 0, N, = 0, f = 1, we have 
N, = z—("%4+%,)(1—f). (42) 
Similarly from (9a) and (1 2m 
N, = ¥,(1—f). (43) 
From (42) and (43), = N,4+N,4+-Y(1-f) (44) 
and from (10), (14), and (42), 
bile ONS F—(H +H MI—f)] 48) 
[Ki —C/3+(C/3)$(f)— Cnd(f)+ Cu, + ¥a)(1—-F ] 
Hence 
(K,—C/8)($'(f)+h+%)+ 
L dp, _ £CR(GEWE N+P gay 








AC df {K, —C/8+-(C/8)d(f)—Cnd(f)+ 
+ Cn(%+%3)(1—f)}* 





"vr = ™ 
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For the form z= (I—f)(1+ @) 
1 kA _ (Ky —€/8)(—14+0+¥, +'¥y)+(C/8)0(%, +) (46b) 
AC| df st {K,—Cn+Cn(4+)}? 


For a tubular charge |@ = 0], when all the charge is burnt, dp,/df has 
same sign as —1+‘"4,+‘. Thus if 


w+, > 1, (47a) 
the maximum pressure occurs before all the charge is burnt, and if 
H+ <1, (47b) 


it occurs at all-burnt stage. 
However, from (42), for a tubular charge 
N, = (—f)—*,—")). 
Thus if our equations are to remain valid (47 b) must be satisfied and for a 
tubular charge in a R.C.L. H./L. gun, the maximum pressure would 
definitely have to occur at the all-burnt stage. 

If (47 a) were satisfied, it would imply, using (40), that the rate of escape 
of the gases from the first chamber is greater than the rate of production 
of gases, which is physically impossible. From (41), ‘4,+‘¥, is proportional 
to S,+<S,, so that in this case the nozzle throat areas would be too big to 
allow any accumulation of gases in the first chamber. 

For an ordinary H./L. gun ‘Y, = 0 and ¥, is of the order of 0-5 and (47 b) 
is satisfied. Thus there also the maximum pressure in the first chamber 
occurs at the all-burnt stage. 

For a cord charge (@ = 1), dp,/df is essentially positive at all-burnt and 
therefore maximum pressure occurs definitely before all-burnt. This result 
is also the same as for orthodox and for ordinary H./L. gun. 

For the quadratic form-function, the condition that maximum pressure 
occurs before the all-burnt stage is 


(K,—C/3)(1—j—¥y) 
(K,—C/8)+-(C/8)(%,+%) 
For a general form-function, the corresponding condition, from (46a), is 


K,| YAW |K— 5 +5 4%) $'(0) > 0. (49) 





@> (48) 


If ¥, = 0, ‘¥, = Y, (48) and (49) reduce to the corresponding conditions 
for ordinary H./L. guns (Kapur (10)). If (49) is satisfied, the maximum 
pressure occurs when 


[Ai —F] + %MI+ SO + HIGH A—s'(N] = 0. (60) 
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The root between 0 and | determines the value of f at which maximum 
pressure occurs in the first chamber and then, substituting this value in (45), 
we get the value of the maximum pressure. 


4.2. Second chamber. Fundamental differential equation 


From (8) and (9), oe AD? (Ps 
ls | = weal p) 











df 
Substituting for p,, p, from (15), (43), and (45), we get 
“| o(f)—(—f (H+) dx) 
df \K,—C/8—C(n—1/3)6(f )+ Cn(¥,+ ¥\(1—f) af} 
_ [ADP 1 ¥,(1—f) 
- in CA K,+Ax—( | P 








. HOR Oe UMS Oa + WON) ni 
$f )—U—f (41+ %2) 
For a tubular charge ¢(f) = 1—f = z, this becomes 


x 





d dX 1+ bz - 
hE id (52) 
dz|1+bz dz X—vz 
where 
» — 1/8—n1—¥,— 
K,—C/6 
pX = K,+Axr (53) 
py = Cy, 


2 A*D? (K,—C/d)*¥, 
~ p CAW,—%—¥)? 
(52) is the same as the corresponding equation of Correr (4) for H./ L. gun, 
only fA and v have to be adjusted. Corner’s series solution for (52) as well 
as his numerical solution in powers of b and v can be used here. 

For a cord charge ¢(f) = 1—f? and (51) can be integrated in series. 
Solutions given by Aggarwal (12) as modified by Kapur (14) can be easily 
adapted for this case. 

Once we know X (or 2) and dX /dz (or dz/df ) as functions of f, we can find 
N, from (42), N, from (43), p, from (45), p, from (15), and 


ic]. 
~ df di ~D™ df 
from the knowledge of p, and dx/df. Also, from (9) and (45), 
1 


= { [K,—C/3+ C/84(f)—Cnd(f)+ Cn +%)(1—f)] 
D CN df) —(¥,4+¥U—f)] 








df. (54) 
f 





), 


) 


) 
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When ¢(f) = 1—f or (1—f)(1+6f), we easily get explicit relations 
between f and f, on evaluating the integral on the right-hand side of (54). 


5. R.C.L. H./L. gun. Non-isothermal model after charge is burnt 
While the isothermal model gives a more or less satisfactory solution 
before the charge is burnt, it breaks down completely after the charge is 
burnt as gases cool off considerably due to expansion. 
Neglecting co-volume terms, the system of equation is given by (6), (8), 
(19), and the following: 


p, K, = CN, RT, (55) 
p{K,+Az] = CN, RT;, (56) 
aN, Ga = |*e p+ Seen] (RT (57) 


aN, 


“A. (58) 


diy _A 
(MT) = — Pp —- S45 yt, 
From (19) and (57), we get on integration 
N, T, 
a= —!. (59) 
(x) Tp 


Again from (6) and (19), we get an integration 

















N,—I,3 = — si5," _—M, si. (60) 
ys 
From (19), (55), and (59), 
aT, 
— = —mT',, (6)) 
bb(S. ~ 
where m ~ See I 
1 
i -2 
Integrating, T, = Tip += , (62) 
where 0, = re 
1B 
From (55), (5 ), and (60), 
,. B 63) 
PB 0; ; ( 





N —tp,)-20-» 
: J (64) 


—-— == | 1+ 
Np 0; 








110 J. N. KAPUR 
From (8) and (58), we get an integration 


; } dx dx\? 
N, T;—N3,2 72.3 = a l(a) ~ (i), | as 
l 


S. _ 
- PprPuel hes 0, (1 +- _ 1\. (65) 


J 


If we use the isothermal model before the charge is burnt, (42), (43) give 
N,. B = 1—Y, —Y, ’ 4 2,.B = Y; . (66) 
Also abies RT, ,-and RT, , by A (65) can be written as 


aile-7a) = ella) —(a),] + 


Aip, pbS, 0 t—t,\*"'7-) 
4M Paw P5201 |() 4 tn _1\ 
2 ©, 
where 7’, is the ratio of the gas temperature to the mean value it had during 
burning. 


—_ 5") 2yy-)) 


(67) 


When ¥, = 0, ‘¥, = V we expect (67) to reduce to an energy equation 
after the charge is burnt for a H./L. gun. We find that the corresponding 
equation of Corner (1, 3) differs from (67) inasmuch as it does not contain 
the last term on the right-hand side of (67), which shows that he has 
neglected there the contribution of elements of gas entering from the first 
chamber into the second. 


Concluding remarks 

We have in the present paper, written a very general system of equations 
of internal ballistics for the recoilless high-low pressure gun and deduced 
the equations for R.C.L., H./L., orthodox guns, and solid-fuel rockets as 
particular cases. Before the all-burnt stage, the isothermal model has been 
studied completely for the R.C.L. H./L. gun, and the solution for this model 
for other guns and rockets can be easily deduced. After the all-burnt stage 
the model breaks down and the non-isothermal model has been studied. 

The isothermal model has been used in Kapur (16) to study the equiva- 
lent non-leaking ballistic problem. The integration of the non-isothermal 
model has also been studied there for the particular case of a tubular charge 
for the linear law and certain interesting results regarding the variation of 
T,, N,, N,T, have been obtained. Once we integrate the non-isothermal 
model for the general law of burning and for the general form-function, 
the main step towards a unified theory would have been taken. The difficul- 
ties are however obvious, as satisfactory solutions, even in particular cases, 
have yet to be discovered. 

Other problems which are under investigation are (i) the effect of 
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Lagrange’s corrections, and (ii) the effect of bore resistance in the present 
unified theory. The second problem is being attacked on the lines used by 
Corner (17) for the orthodox gun. 
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SUMMARY 
Variational expressions for the resonant frequencies of a cavity containing a 
material with tensor electromagnetic properties are derived. The relationship of the 
variational expressions to the Bethe-Schwinger perturbation formula is shown. 


The advantage of this method over the Bethe-Schwinger perturbation method is 
discussed. 


1. Introduction 

MICROWAVE measurements of the electromagnetic properties of materials 
usually involve the effect of a perturbation of a microwave cavity or wave- 
guide by a small sample of the material. Needless to say, the assumption 
that the material possesses a tensor electric permittivity and (or) magnetic 
permeability further complicates an already difficult problem. With the 
increased interest in recent years in the properties and application of ferrites 
in microwave work, one is confronted with the problem of obtaining at least 
approximate solutions to the problem of a cavity or waveguide perturbed 
by an anisotropic obstacle. Perturbation formulae (1) into which one sub- 
stitutes the approximation that the fields equal the unperturbed solution 
or some other suitably chosen fields, have been partially successful in 
obtaining first order results. Perturbation theory however, is not com- 
pletely satisfactory since higher order results are in general tedious to 
obtain, and first order results alone may not be sufficiently accurate. 

This paper is the second (2) of a group of papers in which we concern 
ourselves with the development and application of a general method for 
obtaining approximate solutions to the problems of waveguides and cavities 
containing materials with tensor electromagnetic properties. The method 
is an extension of the work of Julian Schwinger (3) on the theory of obstacles 
in waveguides and cavities to anisotropic media. In each case with the 
introduction of an appropriate dyadic Green’s function we are able to 
obtain a formal solution to the problem in terms of integrals involving the 
field vectors in the perturbing medium. While the resulting integral 
equations are by no means easier to solve, there exists the advantage of 
being able to construct stationary expressions for the quantities of interest 


+ The research reported in this document was supported jointly by the Army, Navy, and 
Air Foree under contract with Massachusetts Institute of Technology. 
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from which the integral equations for the fields within the perturbing 
medium are derivable. Consequently, we have a very powerful method for 
obtaining approximate solutions for these quantities. 

In this paper, for example, we concern ourselves with the problem of a 
cavity perturbed by an anisotropic material. We show that in this case the 
first order variation of (w2—w*), the difference between the square of the 
resonant frequencies of the empty and loaded cavity, is zero with respect 
to similar variations of the fields in the perturbing medium. Berk (4) has 
treated this problem, obtaining variational principles which yield the 
differential equations satisfied by the electromagnetic field. His expres- 
sions, however, are not applicable to dissipative media, and furthermore 
require a knowledge of the fields over the entire cavity or guide. 

While the full use of the variational method consists of the substitution 
of trial functions with unknown variational parameters followed by the 
calculation of the stationary quantity, we may at times obtain good first- 
order results by simply substituting completely determined trial functions. 
In such cases we show how the variational expression reduces to the Bethe— 
Schwinger perturbation formula (5) with the additional advantage of 
yielding the best amplitude for the trial field. 


2. Integral equations and first variational principle for the electro- 
magnetic field 
We consider the problem of an anisotropic obstacle within a cavity free 
of charges and currents. The fields within the cavity satisfy the differential 
equations 


VAE = —jwpn, H—jwp’.H, 
VAH = jwe, E+jwe’.E (1) 
where e, and yu, are respectively the unperturbed constant scalar electric 
permittivity and magnetic permeability. Their choice depends on the 
problem in hand. In the case of obstacles occupying but part of the cavity 
they will be chosen equal to the electric permittivity and magnetic per- 
meability of the empty cavity. The tensor electric permittivity, €, and the 
tensor magnetic permeability, p, of the obstacle are given by 
e=e«I+e’, 
; B= pol+p’, 

where I is the unit dyadic. 

Here e’ and p’ are functions of position, and are zero in the region out- 
side the obstacle. We assume an ei time dependence and observe that 
the fields obey the vector wave equations 


VAVA E—wae, 4, E = wp, €’. E—jw¥ A (p’.H), (2a) 
VAVAH—oe'e,4,H = we, p’.H+jwV A (e’.E). (2b) 


5092.41 I 
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In setting up the variational expression from which one obtains the integral 
equations satisfied by the fields, we shall also require a knowledge of the 
adjoint fields, E' and H', which we choose as the solutions of the equations 


VAE! = jopoH'+joH". p’ 
VAH! = —jwe, E'—jwE*.e’ (3) 
y.B = V.D’ = 0. 


In those problems where € and p are hermitian, we find that H' = H* and 
E’ = E*. In the case of dissipative ferrites, we find 


H'(z, y,z,w) = H(z, y,z, —w) 
and El (z, y,z,w) = E(z, y,z, —w) 


as the off diagonal elements of the antisymmetric magnetic permeability 
tensor are odd in w, and the diagonal elements are even in w (6). 

To solve equation (2a) we introduce the electric dyadic Green’s function 
(3, 7) which satisfies the differential equation 


VAVAN' PN! = [(r—r’) (4) 
and the boundary condition 
nA N'(r{|r’) = 0 


when r lies on the boundary of the cavity. 

Multiplying equation (2a) by N}, = Nta, where a is an arbitrary vector 
introduced in order to simplify the handling of the dyadic Green’s function, 
equation (4) by E(r) and subtracting the result, we find, after integrating 
over the volume of the cavity, that 


E(r’) = wu, [ (e’.E).N*(r |r’) dr—jw | (w’.H).[VAN(r | r’)] dr. 
J (5) 


In those problems where E(r’) is everywhere divergenceless, we shall also 
choose the electric dyadic Green’s function to be divergenceless. 

The normal mode expansion of the dyadic Green’s function as obtained 
by Schwinger (3), is 


i” EX(r)E,(r’) 1 wo fire, (r’) 
Nirir)= > eat BD °) 








n 


where E,, and f,, are respectively the divergenceless and curl-less vector 
eigenfunctions of the empty cavity which satisfy the differential equation 


VAVAA,—RA, = 0 





‘al 
he 


ns 


ty 


on 


or 
on. 





ON THE THEORY OF ANISOTROPIC OBSTACLES IN CAVITIES 115 


and have zero tangential components on the surface of the cavity; A, and 
(2, are normalization constants, k2 the eigenvalues, and 
Aj}(r,w,) = A,(r, —o,,). 

In those regions where E(r) is divergenceless we shall omit the second 
sum in the normal mode expansion of the dyadic Green’s function, except 
for the curl-less eigenfunction, f,, which is also divergenceless. 

To solve for the magnetic field, we may utilize equation (1) or repeat the 
above steps utilizing equation (2b) and the magnetic dyadic Green’s 
function (3). It is readily shown that the two methods are equivalent (3). 

In terms of the normal mode expansion of the dyadic Green’s function 
we have 





; 2. J¢E, (r’) wg 2. ( { fh. e’. ee), - 
- pA bens. o_O 7 
=@ n=0 
where J mn i [w { E}.e’.Edr+w, { H!.p’.H dz| (8) 
and jo, po WH}, = VA E},. 
Taking the curl of equation (7) yields 
p’.A(r’) w, Js H,(r’) _H(r’). (9) 


»  £,» -B) 
Let us note again that in those ani where E(r) is divergenceless we shall 
omit the last sum in equation (7) except for the divergenceless term. 
Proceeding as above we find that 
SJE (r) wg & (JE fede) 0 gy 
(a e) oF aay 


n= 


E'(r) = 





where Je = tela ft e’.E,, dr+w,, | H'.p’.H,, dz| 


H(r).p) Sw, JgtHY(r) 


and = 
i % 2,0 (BP) 


—H'(r). (11) 


From the wave equation satisfied by the electromagnetic fields we find that 
V-[EA (WA E))—E), A (0A E)—jwE) A (p’.H)]+ (kK, —)E,.E 
== wp, wE} .e’. E+, Hj. p’.H). 
It follows therefore that 
Je, = (k2—k?) | Et .E dz/A?. 
We can fix the amplitude of the field within the obstacle by restricting 
ourselves to fields satisfying the condition 


| E.E dr = A, 
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where E, is the unperturbed solution whose perturbation we wish to com- 
t . TI 8 2 © 
pu e 1us Jé ee (k2 —k?), (12) 


and multiplying equation (7) by E’.e’, equation (9) by H'.p’ and adding 
yields Jet A2 = Dé wo, (13) 
where 
De = [ E'.e'.Edr+ [ H'.y’.Hdr+— | Ht. p’.p’.Hdr— 
’ . Mo 
—<h- > Ase ss 2) + 


= Ho 20 


y wy + , + U 
+) cans | | t-« Ear|| [x et, ar|. (14) 


Utilizing equations (10) and (11) it similarly follows that 
J6 AZ = Dh wu. 

Integral equations (7), (9), (10), and (11) for the electromagnetic field are 
by no means easier to solve than the original differential equations. Fortu- 
nately, however, we are able to construct a variational principle for the 
quantity (k2—k?), thus being able to obtain approximate solutions for the 
shift in the resonant frequency. 

To construct such a variational principle, we combine equations (12) 
and (13) to obtain the stationary expression 
Az J¢J¢ 
wg De 
It is readily shown that the first order variation of the right-hand side of 
equation (15) with respect to the fields, subject to the condition that 

Je = Jot = Do w*p,/A3 (16) 

yields the integral equations for the electromagnetic fields within the 
obstacle. Since expression (15) is amplitude-independent this condition is 
automatically satisfied. 


(k3—k?) == (15) 


3. Second variational principle 
A second variational principle may be obtained from the integral equa- 
tions satisfied by the electromagnetic field, which utilize the magnetic 

dyadic Green’s function. Quite analogous to the foregoing we find that 
e’ .E(r’) w, JME, (r’) 


6 4,0 (-P)~ 





E(r’), (17) 


and H(r’) = 








2 mH (r’) ote 5 ( | e1.u-Hdz}e,(') es 


2, a-) 7” 


n=0 





2) 
ug 
3) 


u - 
he 
he 
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where J™ = elo | Et. e’. Edr+w { Hy .p'.H dr], (19) 


and 42 and 7? are the normalization constants of H,(r) and g,,(r) respec- 
tively. The functions g,(r) are the curl-less magnetic eigenfunctions of 
the cavity. In those regions where H is divergenceless we shall omit the 
last sum in equation (18), except for the g, term which is also divergenceless. 
Restricting ourselves to trial functions satisfying the condition 


| H}.H dr = 2, 
we find that (k2—k?) = JR, (20) 
and ITN = JT 'R = wre, D2, (21) 


where Ja" ox “at 7 | [ E. e’.E,, dr+w { Ht. p’-H, dr] 


and 


Dm = [ tart [Be Bart | HY.w Ht ae— 


£9 
MIT IN Jaa w oo + 
“ae ey t Dee ( { ,.u’.-Hdr)( [ Ht.p’.g, dr). 
The variational expression from which equations (17) and (18), and the 


corresponding adjoint equations are derivable is 


Ie Te" 


we, Dm 


(8B) = 





(22) 


We note that had we restricted ourselves to trial fields satisfying the 


condition fz .Edr=A3 (23) 
which we used for the first variational principle, we would have found 
Je = (kW) [ Hy-Hdr)X3 (24) 
So Jt Jo — | t 2 
and Dy = we (ke —k?) = J” f H}.H dr/w Eo- (25) 


Variational expression (22), being amplitude independent, remains 
unaltered. 

The fields which make expression (22) an extremum subject to conditions 
(23) or (25), are the same fields which make expression (15) an extremum. 
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We can therefore combine equations (12) and (24), utilizing equations (8) 
and (19), to find that 


(kk?) eg AR+ ite | H}.H dz| 
= wep Ho(w+)| ( E}.¢’.E dr-+ ( H;}.p’.H dr}, 


l BY.e’.Edr+ ( Ht.p’.Hd 
(w.—w) = w Ia dats or Se (26) 


«,E}.E dr+ | 


by H}.H dr 

This is the well known Bethe—Schwinger perturbation formula (5) general- 
ized to tensor media. The full use of the variational method consists of the 
substitution of trial functions with unknown variational parameters with 
respect to which we then make the variational expression an extremum. 
We may nevertheless at times obtain good first order results by simply 
substituting completely determined trial functions into either equation (12) 
or (20) subject to condition (16) or (21) respectively. Since the variational 
principle permits the computation of the amplitude of the field within the 
perturbing medium, we should expect its results to be better than those 
predicted by the perturbation formula. 
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should like to thank Drs. L. Gold, G.S. Heller, B. Lax, and H. J. Zeiger, and 
Mr. K. J. Button for their interest and help. 
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DETERMINATION OF THE OPTIMUM RESPONSE OF 
LINEAR SYSTEMS 


(ZERO-VELOCITY-ERROR AND ZERO-ACCELERATION- 
ERROR SYSTEMS)+ 


By A. W. BABISTER (The University, Glasgow) 
[Received 29 November 1956] 


SUMMARY 
In (1) two response functions L and L, were defined by the equations 


x 


x 
edr and JL, [ ()’ dr, 
; J \dr. 
0 0 

where ¢ is the error at time 7. In (2) it was shown how these response functions could 
be used to obtain the optimum response for zero-displacement-error systems. The 
analysis is here applied to zero-velocity-error and zero-acceleration-error systems. 
The systems having the optimum response in velocity or in acceleration have a non- 
zero displacement lag. It is shown how the method can be extended to allow for the 
requirement of a zero-displacement-error in the final steady state. 


1. Introduction 
Ix (1) two response functions L and L, were defined by the equations 


L = | edr (1) 


0 


and L,= [ (Zz) a. (2) 
a T 
a 


where ¢ is the error at time r. Formulae were given for L and L, in terms 
of (i) the coefficients of the characteristic equation, and (ii) the frequency 
response spectrum. In (2) it was shown how the response functions could 
be used to obtain the optimum response for zero-displacement-error systems, 
We shall show how the method can be extended to zero-velocity-error and 
zero-acceleration-error systems. 

For simplicity we shall confine our attention to linear systems with one 
degree of freedom. As in (1) and (2), the system considered is one for which 
the equation of motion is 


d"x d"-ly d"-*x 


dx 
a, dr" “T a,-1 dr” 1 Ta,-9 dr"-2 ee + Oy > A oe = f(r), (3) 


* This paper forms part of a thesis for which the degree of Ph.D. of the University of 
Glasgow was awarded. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 1, 1958] 
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where dp, @, @p,..., a, are constants and a, ~ 0. We shall consider two 
cases: (i) where f(r) is a linear function of 7, and (ii) where f(r) is a quadratic 
function. We see that, for a stable system, when the transient motion has 
died away these two cases correspond respectively to systems with (i) con- 
stant velocity, and (ii) constant acceleration. Particular forms of f(r) are 
chosen so that the system ultimately has no position error when subjected 
to either a constant velocity or a constant acceleration input. 


2. Response of zero-velocity-error systems to a linearly increasing 
disturbance (constant velocity input) 


In this case f=0 for7r <0, 


f=fitfor for7r > 0, 
where f, and f, are constants. For simplicity we shall take 
fo _ fi =a, say. 
a @ 


For a stable system when the transient has died away, the steady motion 


is given by oa ae 


and Dz = a, 
where D = d/dr. 

We see that the system may be considered as being subjected to a 
constant velocity input, there being ultimately no position error (since 
x = ar). Such a system is sometimes called a ‘zero-velocity-error system’ 
(see (3), (4), and (5)). 

As in (2) we shall find it convenient to normalize the equation of motion 
(3), but in a different manner from that for the zero-displacement-error 
system. 

Let w, be a frequency defined by 


a, = wi-'*a,, (6) 


so that w, would be the undamped natural frequency of a system with 
ay = 0. We define a new time scale by the relation 


U, = wT. (7) 
From (3)-(7), the normalized equation of motion is 
d"x a*-lz q@*-*% dz. dz 
LP. 5 ——— +9, _» ———_ +... +e —<+-—— +fox 
dun * ? dae t ” 2p * dui du, ? 


=A (141%) (7 > 0), 
8 | 


(m = 0 to n). 





THE OPTIMUM RESPONSE OF LINEAR SYSTEMS 121 


Now the magnitude of the response z is proportional to f,, but the nature 
of the response (percentage overshoot, damping time) is independent of f,. 
We shall take ss (10) 


Then from (5), a=1 and fy, = 4. 
Thus we shall determine the response of the system given by (8) to a 
disturbance f(r) = a,+agr. (11) 


As shown in (1) the response of the system is identical in form with that 
in the free motion with Dz, = —1 and 2, D*x», D®x9,..., D"~!x_ zero where 
the suffix 0 denotes the values at time + = 0. 

For a zero-velocity-error system subjected to a constant velocity input 
we are usually interested (in servomechanism theory) in the velocity 
response at a given time. We shall therefore consider the velocity response 
in the equivalent free motion, finding expressions for 


2 @ 
a dr =| oar (12) 
0 


and (13) 
where = —. (14) 


Second-order zero-velocity-error system 


From (1), with Dz, = —1, 2 = 0, the response functions are given by 


2L,a, = ay, 


(15) 
and 


or in terms of the normalized coefficient r, given by (9), 


2b, = 1+f5. (16) 
wy 
We see that for a system with a given value of w, (i.e. for a second-order 
system, with a given damping factor), L, is a constant and L, is a function 
of r, (which for a mechanical system is proportional to the undamped 
natural frequency). L, is increased by increasing ry. Now for stability 
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ry > 0. Consider the case r, = 0. Equation (8) becomes, using (14), 
dv 
-+v= ] u, > 0). 
du, — (4 ) 
This first-order differential equation shows that v tends monotonically to 
unity as wu, > 00 as shown in Fig. 1 (a). Thus the system r, = 0 is admissible 
and is the optimum system if we are only interested in the response in 


wn 


r (Q) ; r (b) 


lag 


° 
° 


° 
wn 
° 
Ww 


Velocity response 


Displacement 








% 3 i0 0 


Non dimensional time 2, Non dimensional time x, 








Second order zero velocity error system 
Fic. 1. 


velocity, i.e. it is the system for which J, is least for a given value of L,. 
For such a system the displacement x tends to (u,—1)/a), i.e. 

w,(r—2z) > 1, 
giving a displacement lag in the following of a velocity input as shown in 
Fig. 1 (5). 

In a given system it may be either physically impossible or undesirable 
to have r, equal to zero. We have seen that when r, is zero there is a con- 
stant position error (at large times) between the output and the input; this 
is often unacceptable. As shown above, as r, increases L, remains constant 
but ZL, inereases. As with the analysis of zero-displacement-error systems 
(2) we shall choose a value of r, (= }) so that the characteristic equation 
corresponding to (8) has equal roots. The velocity response and displace- 
ment lag are shown in Fig. | plotted against the non-dimensional time w,. 
The velocity response has an overshoot of 14 per cent when r, = }. The 
maximum displacement lag w,(7—.x) is 0-74, the lag becoming less than 0-10 
when u, > 9. Systems having smaller (non-zero) values of r, would have 
a smaller overshoot in velocity but a higher maximum displacement lag. 
Considering the velocity response we are therefore led to the criterion 


0 < rg < 0°25 
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for satisfactory performance for a second-order zero-velocity-error system, 
the precise choice of r, being governed by the maximum acceptable displace- 
ment lag at any time. 


Higher-order zero-velocity-error systems 
The above analysis can easily be extended to higher order systems. 
From (1), with Dx, = —1, we have 


z= Ds, = Px, = ... = “2, = 0, 


21,4 | a @ ...|= Ay A, a, 
Q, 4, + ed 0 ds, a, 
Bese i 0 dy a 


eee eee bie ee OF 


ds pee « |] a, a, 

a, 4; 
a, ‘ ° ne Ay 
Ay 4, an oe | O a, 


or in terms of the normalized coefficients ro, r.,..., 7, given above, 


2 L, W, l Ts rs ', ° e ° = un "s r, e ° ° (17) 
if. To % Mee 2 oi T, 1, 
ee eee | —To YT; 1s 
j | 
ag Me Oa eS fT. 1, 


where the determinants in (17) and (18) are all of order n—1 and r, has 
been put equal to 1. For stable systems of order n (> 2), i.e. systems 
having all the test functions positive, it can be shown that both Z, and 
L, decrease as r, tends to zero. The system with r, = 0 is equivalent to 
an (n—1)th order system in v (= Dax). The velocity response for systems 
making L, a minimum and for systems having equal damping in the least 
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damped modes (referred to later as L,,,;, and L,, systems) follows imme- 
diately from the analysis of (2) (e.g. the L,,;, and I, systems shown in 
Tables 2, 3, and 4 and Fig. 2 of (2)). The displacement lag w,(r—) tends 
to r, (# 0). 

As with the second-order zero-velocity-error system it may be either 
physically impossible or undesirable to have r, equal to zero. Now from 
the analysis of (2) we see that for systems with fzero the L,, system is 
the most satisfactory. To derive a satisfactory velocity-error system with 
ry positive (denoted by L,, system) we replace the zero root of the charac- 
teristic equation of the L,, system by a real root having the same damping 
as the two most lightly damped modes of the L,, system, and normalize 
the resulting equation (with r, = 1). Thus the characteristic equation for 
a third-order zero-velocity-error L,, system is 


AB 2A2--A = (A+1)9A = 0. 


We replace the factor A by A+-1 and normalize the resulting equation; the 
characteristic equation is then 


\8+1-73A2+A+0-19 = (A+0-58)3 = 0 
and the system has three equally damped modes of motion. The values 


of the normalized coefficients and the corresponding roots of the charac- 
teristic equation for higher-order systems are given in Tables | and 2. 


TABLE | 
Values of the normalized coefficients r,, for zero-velocity-error systems 
of order n having three modes with equal damping (L,, systems) 
n To " Ts rs 
73 «I "s 
189 1°53 =I 's 
2°14 2°97 1°52 I " 
3°55 290 366 1°24 I 
TABLE 2 
Roots of the characteristic equation for a zero-velocity-error L,, system 
of order n 
Roots 
0°58, 058, —o58 
—o38, —o°38, —o-38-+1-o1t 
—0°30, —0°'3040°471, —0-304: 1°35! 
0°37, —O17, —O17 10821, —o1741°53i 


As the order of the system increases, the damping of the least-damped 
mode decreases, the percentage overshoot in velocity increases as does the 
maximum displacement lag. 

It is also of interest to normalize the coefficients of the characteristic 
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equation to make the numerical coefficient unity, as with the zero-displace- 
ment-error system discussed in (2). This enables a comparison to be made 
with the results of Whiteley in (3). 


TABLE 3 
Values of the normalized coefficients q,, for zero-velocity-error systems of 
order n having three modes with equal damping (L,, systems) 


n o> 2 ee. 

o> Ree 

4 1 38 46 24 

5 1t 40 Go 59 

6 s 62 108 8&7 7: 18 

Whiteley’s corresponding coefficients are much larger than those shown 

in Table 3, giving modes with very small damping and hence a velocity 
response with a long ‘tail’. This is closely connected with the fact that 
Whiteley’s calculations are based on 10 per cent maximum velocity over- 
shoot whereas in the above analysis much better overall damping has been 
achieved at the expense of a larger overshoot. No curves of displacement 
lags are given in (3) but these would be correspondingly large. Comparing 
the above results with the ITAE curves (5) we find that the velocity 
response of the ITAE system is rather more oscillatory than the L,, curves 
with a maximum overshoot of about the same order. The Z,, curves are 
not unlike those based on binomial filters and are, in general, thought to 
be a reasonable compromise between Whiteley’s curves and the ITAE ones. 


3. Response of zero-acceleration-error systems to a quadratic type 
disturbance (constant acceleration input) 


In this case f=9 forr < 0; 
f=fethitt+hfor? forr > 0, (19) 
where fy, f;, and f, are constants. For simplicity we shall take 


fo_fih 
a, 4, 


= = 5, say. (20) 


1 


For a stable system when the transient has died away, the steady motion 


is given by 


x = $br? 

and Dex = b. 

We see that the system may be considered as being subjected to a constant 

acceleration input, there being ultimately no position error. Such a system 

is sometimes called a ‘zero-acceleration-error system’ (see (3), (4), and (5)). 
The analysis is carried out in a precisely similar manner to that given 

above. We find it convenient to normalize the equation of motion (3) in 
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such a way that the coefficients of d*x/duj and d"x/du} are unity, where 
and 6, = af*a.,. 
From (3), (19)-(22), the normalized equation of motion is 
d"x d®-ly d"-2y _. Bae. d& 
—> 1 8n-1 any 1 on-2 Gong To 188 Ft 
du} du3— du} dug duz 


f. > 
= Pe al uz) (r > 0), (23) 


dx 4 
du, 


+8, 


-89X 


a 


where = _“m_ (m= 0 ton). (24) 
a, w3-™ 


We take J, = Gy. 


As shown in (1) the response of the system is identical in form with that 
in the free motion with D*x, = —1 and x», Dao, D®xpo,..., D"~1xq zero. 

For a zero-acceleration-error system subjected to a constant acceleration 
input we are usually interested in the acceleration response at a given time. 
We shall therefore consider the acceleration response in the equivalent free 
motion, finding values of 


= 2 


d*x\? 
L, — | (7) dr 
0 


x 


3>-\ 2 
and L,; = | (73) dr. (26) 
0 

By a similar analysis to that given above we find that both L, and L, 
decrease as 8, and s, tend to zero. Considering the case s, = s, = 0 we 
find that the relations for L, and L, become identical with those given 
in (2) for L and L, for a zero-displacement-error system of order n—2, 
Ym and wy being replaced by s,,,, and wy. 

For systems with s, = 8s, = 0, neither the displacement lag nor the 
velocity lag tends to zero as 7» 0. This can be remedied as above by 
replacing the two zero roots of the characteristic equation by two real 
roots having the same damping as the two most lightly damped modes 
of the corresponding L,, system and normalizing the resulting equation 
(with s, = 1). The values of the normalized coefficients and the corre- 
sponding roots of the characteristic equation are given in Tables 4 and 5. 

As the order of the system increases (for n > 4) the damping of the least- 
damped mode decreases, 8, and s, only decrease slightly, and the percentage 
overshoot in acceleration increases. 
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TABLE 4 
Values of the normalized coefficients s,, for zero-acceleration-error systems of 
order n having four modes with equal damping (L., systems) 
So Sy Se S3 

0037 O33 +I #I S 

e028 86027 I 163 =«# Ss 

0026 «60°27 1 72 #159 «tf S¢ 

0026 027 +1 «208 «+267 I6I 1 
Note: the third-order system in the above table has, of course, only three 
equally damped modes. 


TABLE 5 
Roots of the characteristic equation for a zero-acceleration-error Lo, system 
of order n 
Roots 


—033 
—o'41, —O°4I 

—0°32, —o732+084i 
—o27+041l, —0'27+1°201 


TABLE 6 
Values of the normalized coefficients q,, for zero-acceleration-error systems of 


order n having four modes with equal damping (Le, systems) 


" Yo 1 G2 qs 
3 3 3 z % 
4 6 4 i %% 
i a ee 4 
6 56 11's 131 gil 30 1 
Note: the third-order system in the above table has, of course, only three 
equally damped modes. 


In Table 6 the coefficients for both third- and fourth-order systems 
are binomial coefficients; this follows immediately from the condition of 
equal damping. Whiteley’s corresponding coefficients (3) are much larger 
than those in Table 6. As stated above, this leads to modes with much 
smaller damping than those shown in Table 5 (allowing for the different 
time scale). The L,, curves are not unlike those based on the binomial 
filters, which are generally thought to be a suitable compromise between 
overshoot and damping. 
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